Splat.h

Go to the documentation of this file.

/** @file */

#ifndef __solve_Splat_h__
#define __solve_Splat_h__

namespace fe
{
namespace ext
{

inline void zero(t_solve_real &a_value)
{
    a_value = 0.0;
}

inline void zero(t_solve_v2 &a_value)
{
    a_value = t_solve_v2(0.0, 0.0);
}

inline void zero(t_solve_v3 &a_value)
{
    a_value = t_solve_v3(0.0, 0.0, 0.0);
}

/// @brief N dimensional state buffer
template<typename T, unsigned int D>
class Splat : public Counted
{
    public:
        typedef Vector<D,int>               t_index;
        typedef t_solve_real                t_span[D];

        class c_gradient
        {
            public:
                c_gradient(void)
                {
                    m_gradient.resize(D);
                    for(unsigned int d = 0; d < D; d++)
                    {
                        zero(m_gradient[d]);
                    }
                }
                std::vector<T>      m_gradient;
        };
        typedef enum
        {
            e_pivot_none = -1,
            e_pivot_within = 0,
            e_pivot_under,
            e_pivot_over,
            e_pivot_count
        } t_pivot_e;
        typedef t_pivot_e                   t_pivot[D];
        typedef c_gradient                  t_gradient;
        typedef std::vector<t_gradient>     t_gradients;
        typedef std::vector<T>              t_cells;
        typedef std::vector<t_solve_real>   t_mass;
        struct t_bnd
        {
            Real    m_lo;
            Real    m_hi;
        };

        void setAll(t_span &a_span, t_solve_real a_value)
        {
            for(unsigned int d = 0; d < D; d++)
            {
                a_span[d] = a_value;
            }
        }


        Splat(void)
        {
            setAll(m_root, 0.0);
            setAll(m_scale, 1.0);
        }
virtual ~Splat(void)
        {
        }
        T                   sample(const typename Splat<T,D>::t_span &a_location);
        void                splat(const t_solve_real &a_mass, const typename Splat<T,D>::t_span &a_location, const T &a_value);
        void                unitSplat(const t_solve_real &a_mass, const typename Splat<T,D>::t_span &a_location, const T &a_value);
        void                massSplat(const t_solve_real &a_mass, const t_index &a_index, const T &a_value);
        void                gradientSplat(const t_solve_real &a_mass, const t_index &a_index, const T &a_value, const t_span &a_direction);
        void                rSplat(unsigned int a_d, typename Splat<T,D>::t_index &a_index, typename Splat<T,D>::t_index &a_lower, float a_mult, t_bnd *a_bnd, const t_solve_real &a_mass, const T &a_value);

        void                create( const t_index &a_count,
                                    const T &a_init);

        const T             &get(const t_index &a_index) const;
        const T             &direct(const t_index &a_index) const;
        const t_solve_real  &getMass(const t_index &a_index) const;
        void                set(const t_index &a_index, const T &a_value);
        void                setMass(const t_index &a_index, const t_solve_real &a_value);

        void                spanize(t_span &a_span, const t_index &a_index);
        t_solve_real        magnitude(const t_span &a_span);

        void                blur(t_solve_real a_factor);
        void                fill(void);
        void                oldstrut(void);
        void                strut(void);
        void                gradient(void);
        void                fillSolve(unsigned int a_aspect);
        void                resplat(Splat<T,D> *a_other_splat);
        void                pull(Splat<T,D> *a_other_splat);
        void                push(Splat<T,D> *a_other_splat);
        void                clone(Splat<T,D> *a_other_splat);

        void                setSpace(   const t_span &a_lo, const t_span &a_hi);
        class   Functor
        {
            public:
                virtual ~Functor(void) {}
                virtual void process(Splat<T,D> &a_splat, t_index &a_index)
                {
                }
        };
        void iterate(Functor &a_functor);

        void                dump(const char *a_filename);
        void                dbsave(const char *a_filename);
        void                dbload(const char *a_filename);
        void                location(   t_span &a_location,
                                        const t_index &a_index) const;

    public:
        t_index         m_count;
        unsigned int    m_totalCount;
        t_cells         m_cells;
        t_mass          m_mass;
        t_span          m_size;
        t_span          m_root;
        t_span          m_scale;
        T               m_init;
        t_gradients     m_gradients;
        t_mass          m_gradient_mass;

        void            aspectFillSolve(unsigned int a_aspect);
        unsigned int    calc_array_index(   const t_index &a_index) const;
        bool            safe(       t_index &a_index) const;
        t_span          &safe(      t_span &a_span) const;
        bool            safe(       unsigned int a_d,
                                    t_index &a_index) const;
        void            unit(       t_span &a_unit_location, const t_span &a_location);
        void            pivot(      t_pivot &a_pivot,
                                    t_index &a_index,
                                    const t_span &a_location) const;
        const t_span    &cellsize(void) const { return m_size; }

        const t_index       &count(void) const { return m_count; }
    private:
        T                   &access(const t_index &a_index);
        T                   rLinearSample(unsigned int a_d, typename Splat<T,D>::t_index &a_index, typename Splat<T,D>::t_index &a_lower, typename Splat<T,D>::t_pivot &a_pivot, float a_mult, t_bnd *a_bnd);
        void rIterate(unsigned int a_d, t_index &a_index, Functor &a_functor);

        unsigned int        totalCount(void) const { return m_totalCount; }
        void                index(      t_index &a_index,
                                        const t_span &a_location) const;
};

template<typename T, unsigned int D>
void Splat<T,D>::splat(const t_solve_real &a_mass, const typename Splat<T,D>::t_span &a_location, const T &a_value)
{
    t_bnd bnd[D];

    Splat<T,D>::t_span unit_location;
    unit(unit_location, a_location);

    t_pivot i_pivot;
    t_index i_lower;
    pivot(i_pivot, i_lower, unit_location);

    for(unsigned int i = 0; i < D; i++)
    {
        if(i_pivot[i] == e_pivot_within)
        {
            bnd[i].m_lo = unit_location[i]/cellsize()[i] - (Real)i_lower[i];
            bnd[i].m_hi = 1.0 - bnd[i].m_lo;
        }
        else if(i_pivot[i] == e_pivot_over)
        {
            //bnd[i].m_hi = 0.0;
            //bnd[i].m_lo = 1.0;
            return;
        }
        else if(i_pivot[i] == e_pivot_under)
        {
            //bnd[i].m_hi = 1.0;
            //bnd[i].m_lo = 0.0;
            return;
        }
    }

    Real r = 1.0;

    t_index index;
    rSplat(D-1, index, i_lower, r, bnd, a_mass, a_value);
}

template<typename T, unsigned int D>
void Splat<T,D>::unitSplat(const t_solve_real &a_mass, const typename Splat<T,D>::t_span &a_location, const T &a_value)
{
    t_bnd bnd[D];

    t_pivot i_pivot;
    t_index i_lower;
    pivot(i_pivot, i_lower, a_location);

    for(unsigned int i = 0; i < D; i++)
    {
        if(i_pivot[i] == e_pivot_within)
        {
            bnd[i].m_lo = a_location[i]/cellsize()[i] - (Real)i_lower[i];
            bnd[i].m_hi = 1.0 - bnd[i].m_lo;
        }
        else if(i_pivot[i] == e_pivot_over)
        {
            bnd[i].m_hi = 0.0;
            bnd[i].m_lo = 1.0;
        }
        else if(i_pivot[i] == e_pivot_under)
        {
            bnd[i].m_hi = 1.0;
            bnd[i].m_lo = 0.0;
        }
    }

    Real r = 1.0;

    t_index index;
    rSplat(D-1, index, i_lower, r, bnd, a_mass, a_value);
}

template<typename T, unsigned int D>
void Splat<T,D>::massSplat(const t_solve_real &a_mass, const t_index &a_index, const T &a_value)
{
    if(a_mass < 1.0e-12) { return; }
    t_index safe_index(a_index);
    safe(safe_index);
    unsigned int i = calc_array_index(safe_index);
    t_solve_real sum = a_mass + m_mass[i];
    m_cells[i] = (m_mass[i] / sum) * m_cells[i] + (a_mass / sum) * a_value;
    m_mass[i] += a_mass;
}

template<typename T, unsigned int D>
inline void Splat<T,D>::location(t_span &a_location, const t_index &a_index) const
{
    for(unsigned int i = 0; i < D; i++)
    {
        a_location[i] = m_root[i] + (Real)a_index[i] * cellsize()[i] / m_scale[i];
    }
}

template<typename T, unsigned int D>
void Splat<T,D>::gradientSplat(const t_solve_real &a_mass, const t_index &a_index, const T &a_value, const t_span &a_unit_dir)
{
    if(a_mass < 1.0e-12) { return; }

    t_index safe_index(a_index);
    safe(safe_index);
    unsigned int i_cell = calc_array_index(safe_index);
    t_gradient &gradient = m_gradients[i_cell];

    t_solve_real sum = a_mass + m_gradient_mass[i_cell];
    t_solve_real sum_ratio = a_mass / sum;
    t_solve_real sum_m_ratio = m_gradient_mass[i_cell] / sum;
    m_gradient_mass[i_cell] += a_mass;

    for(unsigned int d = 0; d < D; d++)
    {
        gradient.m_gradient[d] = sum_m_ratio * gradient.m_gradient[d] + sum_ratio * (a_value * a_unit_dir[d]);
    }
}

template<typename T, unsigned int D>
void Splat<T,D>::setSpace(  const t_span &a_lo, const t_span &a_hi)
{
    for(unsigned int d = 0; d < D; d++)
    {
        m_root[d] = a_lo[d];
        m_scale[d] = 1.0 / (a_hi[d] - a_lo[d]);
    }
}


template<typename T, unsigned int D>
void Splat<T,D>::unit(  t_span &a_unit_location, const t_span &a_location)
{
    for(unsigned int i = 0; i < D; i++)
    {
        a_unit_location[i] = (a_location[i] - m_root[i]) * m_scale[i];
    }
}

template<typename T, unsigned int D>
inline T &Splat<T,D>::access(const t_index &a_index)
{
    return m_cells[calc_array_index(a_index)];
}

template<typename T, unsigned int D>
inline void Splat<T,D>::index(t_index &a_index, const t_span &a_location) const
{
    for(unsigned int i = 0; i < D; i++)
    {
        a_index[i] = (int)(a_location[i] / m_size[i]);
    }
}

template<typename T, unsigned int D>
inline void Splat<T,D>::pivot(t_pivot &a_pivot, t_index &a_index, const t_span &a_location) const
{
    for(unsigned int i = 0; i < D; i++)
    {
        a_index[i] = (int)((a_location[i]) / m_size[i]);

        if((int)a_index[i] < 0)
        {
            a_pivot[i] = e_pivot_under;
            a_index[i] = 0;
        }
        else if((int)a_index[i] >= (int)(m_count[i])-1)
        {
            a_pivot[i] = e_pivot_over;
            a_index[i] = m_count[i]-1;
        }
        else
        {
            a_pivot[i] = e_pivot_within;
        }
    }
}

template<typename T, unsigned int D>
inline void Splat<T,D>::spanize(t_span &a_span, const t_index &a_index)
{
    for(unsigned int i = 0; i < D; i++)
    {
        a_span[i] = (Real)a_index[i] * cellsize()[i];
    }
}

template<typename T, unsigned int D>
inline t_solve_real Splat<T,D>::magnitude(const t_span &a_span)
{
    t_solve_real m = 0.0;
    for(unsigned int d = 0; d < D; d++)
    {
        m += a_span[d]*a_span[d];
    }
    return sqrt(m);
}

template<typename T, unsigned int D>
void Splat<T,D>::rSplat(unsigned int a_d, typename Splat<T,D>::t_index &a_index, typename Splat<T,D>::t_index &a_lower, float a_mult, t_bnd *a_bnd, const t_solve_real &a_mass, const T &a_value)
{
    Real m;
    a_index[a_d] = a_lower[a_d];
    safe(a_index);
    m = a_bnd[a_d].m_hi * a_mult;
    if(m > 0.0)
    {
        if(a_d == 0)
        {
            massSplat(m * a_mass, a_index, a_value);
        }
        else
        {
            rSplat(a_d-1, a_index, a_lower, m, a_bnd, a_mass, a_value);
        }
    }

    a_index[a_d] = a_lower[a_d] + 1;
    safe(a_index);
    m = a_bnd[a_d].m_lo * a_mult;
    if(m > 0.0)
    {
        if(a_d == 0)
        {
            massSplat(m * a_mass, a_index, a_value);
        }
        else
        {
            rSplat(a_d-1, a_index, a_lower, m, a_bnd, a_mass, a_value);
        }
    }
}

inline t_solve_real &value(t_solve_real &a_value, unsigned int a_aspect)
{
    return a_value;
}

inline t_solve_real &value(t_solve_v2 &a_value, unsigned int a_aspect)
{
    return a_value[a_aspect];
}

inline t_solve_real &value(t_solve_v3 &a_value, unsigned int a_aspect)
{
    return a_value[a_aspect];
}

template<typename T, unsigned int D>
T Splat<T,D>::rLinearSample(unsigned int a_d, typename Splat<T,D>::t_index &a_index, typename Splat<T,D>::t_index &a_lower, typename Splat<T,D>::t_pivot &a_pivot, float a_mult, t_bnd *a_bnd)
{
    T accum;
    zero(accum);

    bool ok = false;
    Real m;
    a_index[a_d] = a_lower[a_d];
    ok = safe(a_index);
    m = a_bnd[a_d].m_hi * a_mult;
    //if(m > 0.0)
    {
        if(a_d == 0)
        {
            accum += m * direct(a_index);
        }
        else
        {
            accum += rLinearSample(a_d-1, a_index, a_lower, a_pivot, m, a_bnd);
        }
    }


    a_index[a_d] = a_lower[a_d] + 1;
    ok = safe(a_index);
    m = a_bnd[a_d].m_lo * a_mult;

    //if(m > 0.0)
    {
        if(a_d == 0)
        {
            accum += m * direct(a_index);
        }
        else
        {
            accum += rLinearSample(a_d-1, a_index, a_lower, a_pivot, m, a_bnd);
        }
    }

    return accum;
}

/** sample a Continuum by returning the a linear interpolated value.
    @relates Continuum */
template<typename T, unsigned int D>
T Splat<T,D>::sample(const typename Splat<T,D>::t_span &a_location)
{
    t_bnd bnd[D];

    Splat<T,D>::t_span unit_location;
    unit(unit_location, a_location);

    t_pivot i_pivot;
    t_index i_lower;
    pivot(i_pivot, i_lower, unit_location);

    for(unsigned int i = 0; i < D; i++)
    {
        if(i_pivot[i] == e_pivot_within)
        {
            bnd[i].m_lo = unit_location[i]/cellsize()[i] - (Real)i_lower[i];
            bnd[i].m_hi = 1.0 - bnd[i].m_lo;
        }
        else if(i_pivot[i] == e_pivot_over)
        {
            i_lower[i] -= 1;
            bnd[i].m_hi = 1.0;
            bnd[i].m_lo = 0.0;
        }
        else if(i_pivot[i] == e_pivot_under)
        {
            bnd[i].m_hi = 1.0;
            bnd[i].m_lo = 0.0;
        }
    }

    Real r = 1.0;
    t_index index;
    return rLinearSample(D-1, index, i_lower, i_pivot, r, bnd);
}

inline t_solve_v2 sample(Splat<t_solve_v2,3> &a_splat, const typename Splat<t_solve_v2,3>::t_span &a_location)
{
    Splat<t_solve_v2,3>::t_bnd bnd[3];
    Splat<t_solve_v2,3>::t_span unit_location;
    a_splat.unit(unit_location, a_location);

    Splat<t_solve_v2,3>::t_pivot i_pivot;
    Splat<t_solve_v2,3>::t_index i_lower;
    a_splat.pivot(i_pivot, i_lower, unit_location);

    for(unsigned int i = 0; i < 3; i++)
    {
        if(i_pivot[i] == Splat<t_solve_v2,3>::e_pivot_within)
        {
            bnd[i].m_lo = unit_location[i]/a_splat.cellsize()[i] - (Real)i_lower[i];
            bnd[i].m_hi = 1.0 - bnd[i].m_lo;
        }
        else if(i_pivot[i] == Splat<t_solve_v2,3>::e_pivot_over)
        {
            i_lower[i] -= 1;
            bnd[i].m_hi = 1.0;
            bnd[i].m_lo = 0.0;
        }
        else if(i_pivot[i] == Splat<t_solve_v2,3>::e_pivot_under)
        {
            bnd[i].m_hi = 1.0;
            bnd[i].m_lo = 0.0;
        }
    }

    t_solve_v2 accum(0.0,0.0);
    Splat<t_solve_v2,3>::t_index i_sample;
    i_sample[0]=i_lower[0];i_sample[1]=i_lower[1];i_sample[2]=i_lower[2];
    accum += a_splat.direct(i_sample) * bnd[0].m_hi * bnd[1].m_hi * bnd[2].m_hi;

    i_sample[0]=i_lower[0];i_sample[1]=i_lower[1];i_sample[2]=i_lower[2]+1;
    accum += a_splat.direct(i_sample) * bnd[0].m_hi * bnd[1].m_hi * bnd[2].m_lo;

    i_sample[0]=i_lower[0];i_sample[1]=i_lower[1]+1;i_sample[2]=i_lower[2]+1;
    accum += a_splat.direct(i_sample) * bnd[0].m_hi * bnd[1].m_lo * bnd[2].m_lo;

    i_sample[0]=i_lower[0]+1;i_sample[1]=i_lower[1]+1;i_sample[2]=i_lower[2]+1;
    accum += a_splat.direct(i_sample) * bnd[0].m_lo * bnd[1].m_lo * bnd[2].m_lo;

    i_sample[0]=i_lower[0]+1;i_sample[1]=i_lower[1]+1;i_sample[2]=i_lower[2];
    accum += a_splat.direct(i_sample) * bnd[0].m_lo * bnd[1].m_lo * bnd[2].m_hi;

    i_sample[0]=i_lower[0]+1;i_sample[1]=i_lower[1];i_sample[2]=i_lower[2];
    accum += a_splat.direct(i_sample) * bnd[0].m_lo * bnd[1].m_hi * bnd[2].m_hi;

    i_sample[0]=i_lower[0]+1;i_sample[1]=i_lower[1];i_sample[2]=i_lower[2]+1;
    accum += a_splat.direct(i_sample) * bnd[0].m_lo * bnd[1].m_hi * bnd[2].m_lo;

    i_sample[0]=i_lower[0];i_sample[1]=i_lower[1]+1;i_sample[2]=i_lower[2];
    accum += a_splat.direct(i_sample) * bnd[0].m_hi * bnd[1].m_lo * bnd[2].m_hi;
    return accum;
}

inline t_solve_v2 sample(Splat<t_solve_v2,2> &a_splat, const typename Splat<t_solve_v2,2>::t_span &a_location)
{
    Splat<t_solve_v2,2>::t_bnd bnd[2];
    Splat<t_solve_v2,2>::t_span unit_location;
    a_splat.unit(unit_location, a_location);

    Splat<t_solve_v2,2>::t_pivot i_pivot;
    Splat<t_solve_v2,2>::t_index i_lower;
    a_splat.pivot(i_pivot, i_lower, unit_location);

    for(unsigned int i = 0; i < 2; i++)
    {
        if(i_pivot[i] == Splat<t_solve_v2,2>::e_pivot_within)
        {
            bnd[i].m_hi = unit_location[i]/a_splat.cellsize()[i] - (Real)i_lower[i];
            bnd[i].m_lo = 1.0 - bnd[i].m_hi;
        }
        else if(i_pivot[i] == Splat<t_solve_v2,2>::e_pivot_over)
        {
            i_lower[i] -= 1;
            bnd[i].m_hi = 0.0;
            bnd[i].m_lo = 1.0;
        }
        else if(i_pivot[i] == Splat<t_solve_v2,2>::e_pivot_under)
        {
            bnd[i].m_hi = 0.0;
            bnd[i].m_lo = 1.0;
        }
    }

    t_solve_v2 accum(0.0,0.0);
    Splat<t_solve_v2,2>::t_index i_sample;
    i_sample[0]=i_lower[0];i_sample[1]=i_lower[1];
    accum += a_splat.direct(i_sample) * bnd[0].m_lo * bnd[1].m_lo;

    i_sample[0]=i_lower[0];i_sample[1]=i_lower[1]+1;
    accum += a_splat.direct(i_sample) * bnd[0].m_lo * bnd[1].m_hi;

    i_sample[0]=i_lower[0]+1;i_sample[1]=i_lower[1]+1;
    accum += a_splat.direct(i_sample) * bnd[0].m_hi * bnd[1].m_hi;

    i_sample[0]=i_lower[0]+1;i_sample[1]=i_lower[1];
    accum += a_splat.direct(i_sample) * bnd[0].m_hi * bnd[1].m_lo;

    return accum;
}

template<typename T, unsigned int D>
class SubFill : virtual public Splat<T,D>::Functor
{
    public:
        SubFill(typename Splat<T,D>::t_index &a_index, const T &a_value, t_solve_real a_mass, std::vector<t_solve_v2> *a_data)
        {
            m_index = a_index;
            m_value = a_value;
            m_mass = a_mass;
            m_data = a_data;
        }
        virtual ~SubFill(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            t_solve_real sum = 0;
            //t_solve_real axis_factor = 1.0;
            //unsigned off_cnt = 0;
            for(unsigned int i = 0; i < D; i++)
            {
                t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                sum += a*a;
            }
            t_solve_real distance = sqrt(sum);
            if(distance < (*m_data)[i_cell][0])
            {
                t_solve_real dial = 0.1;
                if(a_splat.m_mass[i_cell] <= dial && m_mass > dial)
                {
                    (*m_data)[i_cell][0] = distance;
                    a_splat.m_cells[i_cell] = m_value;
                    a_splat.m_mass[i_cell] = dial;
                }
            }
        }
        typename Splat<T,D>::t_index m_index;
        T m_value;
        std::vector<t_solve_v2> *m_data;
        t_solve_real m_mass;
};

template<typename T, unsigned int D>
class Fill : virtual public Splat<T,D>::Functor
{
    public:
        Fill(typename Splat<T,D>::t_cells *a_cells, typename Splat<T,D>::t_mass *a_mass)
        {
            m_src_cells = a_cells;
            m_src_mass = a_mass;
            m_data.resize(m_src_mass->size());
            for(unsigned int i = 0; i < m_data.size(); i++)
            {
                m_data[i] = t_solve_v2(1.0e8,1.0e8);
            }
        }
        virtual ~Fill(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            T &src_value = (*m_src_cells)[i_cell];
            t_solve_real src_mass = (*m_src_mass)[i_cell];
            SubFill<T,D> subfill(a_index, src_value, src_mass, &m_data);
            a_splat.iterate(subfill);
        }
        typename Splat<T,D>::t_cells *m_src_cells;
        typename Splat<T,D>::t_mass *m_src_mass;
        std::vector<t_solve_v2> m_data;
};

template<typename T, unsigned int D>
void Splat<T,D>::fill(void)
{
    t_cells src_cells;
    t_mass src_mass;
    src_cells = m_cells;
    src_mass = m_mass;

    Fill<T,D> fill(&src_cells, &src_mass);
    iterate(fill);
}


template<typename T, unsigned int D>
class SubStrut : virtual public Splat<T,D>::Functor
{
    public:
        SubStrut(typename Splat<T,D>::t_index &a_index, const T &a_value, typename Splat<T,D>::t_mass *a_mass, std::vector<t_solve_v2> *a_data, typename Splat<T,D>::t_index *a_closest0, typename Splat<T,D>::t_index *a_closest1)
        {
            m_index = a_index;
            m_value = a_value;
            m_mass = a_mass;
            m_data = a_data;
            m_closest0 = a_closest0;
            m_closest1 = a_closest1;
        }
        virtual ~SubStrut(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            t_solve_real sum = 0;
            t_solve_real closest0 = 0;
            t_solve_real closest1 = 0;
            //t_solve_real axis_factor = 1.0;
            //unsigned off_cnt = 0;
            for(unsigned int i = 0; i < D; i++)
            {
                t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                sum += a*a;

                a = (t_solve_real)((*m_closest0)[i] - m_index[i]) * a_splat.m_size[i];
                closest0 += a*a;

                a = (t_solve_real)((*m_closest1)[i] - m_index[i]) * a_splat.m_size[i];
                closest1 += a*a;
            }


            if((*m_mass)[i_cell] > 1.0e-6)
            {
                if(sum < closest0)
                {
                    (*m_closest1) = (*m_closest0);
                    (*m_closest0) = a_index;

                    sum = 0.0; closest0 = 0.0; closest1 = 0.0;
                    for(unsigned int i = 0; i < D; i++)
                    {
                        t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                        sum += a*a;
                        a = (t_solve_real)((*m_closest0)[i] - m_index[i]) * a_splat.m_size[i];
                        closest0 += a*a;
                        a = (t_solve_real)((*m_closest1)[i] - m_index[i]) * a_splat.m_size[i];
                        closest1 += a*a;
                    }
                }
                else if(sum < closest1)
                {
                    (*m_closest1) = a_index;
                    sum = 0.0; closest0 = 0.0; closest1 = 0.0;
                    for(unsigned int i = 0; i < D; i++)
                    {
                        t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                        sum += a*a;
                        a = (t_solve_real)((*m_closest0)[i] - m_index[i]) * a_splat.m_size[i];
                        closest0 += a*a;
                        a = (t_solve_real)((*m_closest1)[i] - m_index[i]) * a_splat.m_size[i];
                        closest1 += a*a;
                    }
                }
            }
        }
        typename Splat<T,D>::t_index m_index;
        T m_value;
        std::vector<t_solve_v2> *m_data;
        typename Splat<T,D>::t_mass *m_mass;
        typename Splat<T,D>::t_index *m_closest0;
        typename Splat<T,D>::t_index *m_closest1;
};

template<typename T, unsigned int D>
class Closest : virtual public Splat<T,D>::Functor
{
    public:
        Closest(typename Splat<T,D>::t_index &a_index, const T &a_value, typename Splat<T,D>::t_mass *a_mass, std::vector<t_solve_v2> *a_data, typename Splat<T,D>::t_index *a_closest0)
        {
            m_index = a_index;
            m_value = a_value;
            m_mass = a_mass;
            m_data = a_data;
            m_closest0 = a_closest0;
        }
        virtual ~Closest(void)
        {
        }
        void exclude(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            m_exclude.push_back(i_cell);
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            for(unsigned int i = 0; i < m_exclude.size(); i++)
            {
                if(i_cell == m_exclude[i])
                {
                    return;
                }
            }
            t_solve_real sum = 0;
            t_solve_real closest0 = 0;
            for(unsigned int i = 0; i < D; i++)
            {
                t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                sum += a*a;

                a = (t_solve_real)((*m_closest0)[i] - m_index[i]) * a_splat.m_size[i];
                closest0 += a*a;
            }


            if((*m_mass)[i_cell] > 1.0e-6)
            {
                if(sum < closest0)
                {
                    (*m_closest0) = a_index;

                    sum = 0.0; closest0 = 0.0;
                    for(unsigned int i = 0; i < D; i++)
                    {
                        t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                        sum += a*a;
                        a = (t_solve_real)((*m_closest0)[i] - m_index[i]) * a_splat.m_size[i];
                        closest0 += a*a;
                    }
                }
            }
        }
        typename Splat<T,D>::t_index m_index;
        std::vector<unsigned int> m_exclude;
        T m_value;
        std::vector<t_solve_v2> *m_data;
        typename Splat<T,D>::t_mass *m_mass;
        typename Splat<T,D>::t_index *m_closest0;
};

template<typename T, unsigned int D>
class ClosestMate : virtual public Splat<T,D>::Functor
{
    public:
        ClosestMate(typename Splat<T,D>::t_index &a_index, const T &a_value, typename Splat<T,D>::t_mass *a_mass, std::vector<t_solve_v2> *a_data, typename Splat<T,D>::t_index *a_closest, typename Splat<T,D>::t_index *a_closest_mate)
        {
            m_index = a_index;
            m_value = a_value;
            m_mass = a_mass;
            m_data = a_data;
            m_closest = a_closest;
            m_closest_mate = a_closest_mate;
        }
        virtual ~ClosestMate(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);

            typename Splat<T,D>::t_span M;
            a_splat.spanize(M, (*m_closest) - m_index);
            typename Splat<T,D>::t_span C;
            a_splat.spanize(C, a_index - m_index);

            t_solve_real DOT = 0.0;
            t_solve_real MM = 0.0;
            t_solve_real CC = 0.0;
            for(unsigned int d = 0; d < D; d++)
            {
                MM += M[d]*M[d];
                CC += C[d]*C[d];
            }
            MM = sqrt(MM);
            CC = sqrt(CC);
if(MM < 1.0e-8) { return; }
if(CC < 1.0e-8) { return; }

            for(unsigned int d = 0; d < D; d++)
            {
                DOT += M[d]/MM * C[d]/MM;
            }

            if(DOT >= -1.0e-8) { return; }
            //if(DOT >= 0.5) { return; }

            t_solve_real sum = 0;
            t_solve_real closest_mate = 0;
            for(unsigned int i = 0; i < D; i++)
            {
                t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                sum += a*a;

                a = (t_solve_real)((*m_closest_mate)[i] - m_index[i]) * a_splat.m_size[i];
                closest_mate += a*a;
            }


            if((*m_mass)[i_cell] > 1.0e-6)
            {
                if(sum < closest_mate)
                {
                    (*m_closest_mate) = a_index;
                }
            }
        }
        typename Splat<T,D>::t_index m_index;
        T m_value;
        std::vector<t_solve_v2> *m_data;
        typename Splat<T,D>::t_mass *m_mass;
        typename Splat<T,D>::t_index *m_closest;
        typename Splat<T,D>::t_index *m_closest_mate;
};

template<typename T, unsigned int D>
class Strut : virtual public Splat<T,D>::Functor
{
    public:
        Strut(typename Splat<T,D>::t_cells *a_cells, typename Splat<T,D>::t_mass *a_mass)
        {
            m_src_cells = a_cells;
            m_src_mass = a_mass;
            m_data.resize(m_src_mass->size());
            for(unsigned int i = 0; i < m_data.size(); i++)
            {
                m_data[i] = t_solve_v2(1.0e8,1.0e8);
            }
        }
        virtual ~Strut(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            T &src_value = (*m_src_cells)[i_cell];
            t_solve_real src_mass = (*m_src_mass)[i_cell];
#if 1
if(src_mass > 1.0e-8)
{
    a_splat.set(a_index, src_value);
    a_splat.setMass(a_index, 1.0);
    return;
}
#endif

            for(unsigned int d = 0; d < D; d++)
            {
                m_closest0[d] = 1000;
                m_closest1[d] = 1000;
                m_closest2[d] = 1000;
                m_closest0mate[d] = 1000;
                m_closest1mate[d] = 1000;
                m_closest2mate[d] = 1000;
            }
            //SubStrut<T,D> substrut(a_index, src_value, m_src_mass, &m_data, &m_closest0, &m_closest1);
            //a_splat.iterate(substrut);

            {
                Closest<T,D> closest(a_index, src_value, m_src_mass, &m_data, &m_closest0);
                closest.exclude(a_splat, a_index);
                a_splat.iterate(closest);
            }
            {
                ClosestMate<T,D> closestmate(a_index, src_value, m_src_mass, &m_data, &m_closest0, &m_closest0mate);
                a_splat.iterate(closestmate);
            }

            {
                Closest<T,D> closest(a_index, src_value, m_src_mass, &m_data, &m_closest1);
                closest.exclude(a_splat, a_index);
                closest.exclude(a_splat, m_closest0);
                closest.exclude(a_splat, m_closest0mate);
                a_splat.iterate(closest);
            }

            if(m_closest0mate[0] != 1000)
            {
                ClosestMate<T,D> closestmate(a_index, src_value, m_src_mass, &m_data, &m_closest0mate, &m_closest1mate);
                a_splat.iterate(closestmate);
            }






            {
                Closest<T,D> closest(a_index, src_value, m_src_mass, &m_data, &m_closest2);
                closest.exclude(a_splat, a_index);
                closest.exclude(a_splat, m_closest0mate);
                a_splat.iterate(closest);
            }
            if(m_closest2[0] != 1000)
            {
                ClosestMate<T,D> closestmate(a_index, src_value, m_src_mass, &m_data, &m_closest2, &m_closest2mate);
                a_splat.iterate(closestmate);
            }



            t_solve_real closest0 = 0;
            t_solve_real closest1 = 0;
            t_solve_real closest2 = 0;
            t_solve_real closest0mate = 0;
            t_solve_real closest1mate = 0;
            t_solve_real closest2mate = 0;
            for(unsigned int i = 0; i < D; i++)
            {
                t_solve_real a;
                a = (t_solve_real)(m_closest0[i] - a_index[i]) * a_splat.m_size[i];
                closest0 += a*a;

                a = (t_solve_real)(m_closest1[i] - a_index[i]) * a_splat.m_size[i];
                closest1 += a*a;

                a = (t_solve_real)(m_closest2[i] - a_index[i]) * a_splat.m_size[i];
                closest2 += a*a;

                a = (t_solve_real)(m_closest0mate[i] - a_index[i]) * a_splat.m_size[i];
                closest0mate += a*a;

                a = (t_solve_real)(m_closest1mate[i] - a_index[i]) * a_splat.m_size[i];
                closest1mate += a*a;

                a = (t_solve_real)(m_closest2mate[i] - a_index[i]) * a_splat.m_size[i];
                closest2mate += a*a;
            }

            typename Splat<T,D>::t_span U;
            a_splat.spanize(U, m_closest1 - m_closest0);
            typename Splat<T,D>::t_span P;
            a_splat.spanize(P, a_index - m_closest0);

            //typename Splat<T,D>::t_span C0;
            //a_splat.spanize(C0, m_closest0);
            //typename Splat<T,D>::t_span C1;
            //a_splat.spanize(C1, m_closest1);
            typename Splat<T,D>::t_span AI;
            a_splat.spanize(AI, a_index);

            t_solve_real c0z = 0.0;
            t_solve_real L = 0.0;
            for(unsigned int d = 0; d < D; d++)
            {
                c0z += P[d]*P[d];
                L += U[d] * U[d];
            }
            c0z = sqrt(c0z);
            L = sqrt(L);

            L = 1.0 / L;

#if 0
            t_solve_real DOT = 0.0;
            for(unsigned int d = 0; d < D; d++)
            {
                DOT += P[d] * U[d];
            }
            DOT *= L*L;
#endif

            T c0 = (*m_src_cells)[a_splat.calc_array_index(m_closest0)];
            if(m_closest0mate[0] == 1000
                || m_closest1mate[0] == 1000 )
            {
                T G;
                zero(G);
                for(unsigned int d = 0; d < D; d++)
                {
                    G = G + a_splat.m_gradients[a_splat.calc_array_index(m_closest0)].m_gradient[d] * P[d];
                }
                a_splat.set(a_index, c0 + G);
                //a_splat.set(a_index, c0);
                closest0 = sqrt(closest0);
                a_splat.setMass(a_index, 1.0-closest0);
                //a_splat.setMass(a_index, 0.01);
            }
            else
            {
                T c0m = (*m_src_cells)[a_splat.calc_array_index(m_closest0mate)];
                closest0 = sqrt(closest0);
                closest0mate = sqrt(closest0mate);

                T c1 = (*m_src_cells)[a_splat.calc_array_index(m_closest1)];
                T c1m = (*m_src_cells)[a_splat.calc_array_index(m_closest1mate)];
                closest1 = sqrt(closest1);
                closest1mate = sqrt(closest1mate);

                //T c2 = (*m_src_cells)[a_splat.calc_array_index(m_closest2)];
                //T c2m = (*m_src_cells)[a_splat.calc_array_index(m_closest2mate)];
                closest2 = sqrt(closest2);
                closest2mate = sqrt(closest2mate);
#if 0
                T G0; zero(G0);
                T G1; zero(G1);
                for(unsigned int d = 0; d < D; d++)
                {
                    G0 = G0 + a_splat.m_gradients[a_splat.calc_array_index(m_closest0)].m_gradient[d] * P[d];
                    G1 = G1 + a_splat.m_gradients[a_splat.calc_array_index(m_closest1)].m_gradient[d] * P[d];
                }
#endif
                t_solve_real w0 = closest0mate / (closest0 + closest0mate);
                t_solve_real w0m = 1.0 - w0;
                t_solve_real w1 = closest1mate / (closest1 + closest1mate);
                t_solve_real w1m = 1.0 - w1;
                //t_solve_real w2 = closest2mate / (closest2 + closest2mate);
                //t_solve_real w2m = 1.0 - w2;

                t_solve_real W0 = closest1 / (closest0 + closest1);
                t_solve_real W1 = 1.0 - W0;

                T C0 = w0*c0 + w0m*c0m;
                T C1 = w1*c1 + w1m*c1m;

                T C = W0*C0 + W1*C1;
C = C0;
                a_splat.set(a_index, C);
                //a_splat.setMass(a_index, w0);
                a_splat.setMass(a_index, 1.0 - closest0);
//fe_fprintf(stderr, "%g\n", w0);
                //a_splat.set(a_index, c0);
                //a_splat.set(a_index, w0*(c0 + G0) + w1*(c1 + G1));
            }
        }
        typename Splat<T,D>::t_cells *m_src_cells;
        typename Splat<T,D>::t_mass *m_src_mass;
        typename Splat<T,D>::t_index m_closest0;
        typename Splat<T,D>::t_index m_closest1;
        typename Splat<T,D>::t_index m_closest2;
        typename Splat<T,D>::t_index m_closest0mate;
        typename Splat<T,D>::t_index m_closest1mate;
        typename Splat<T,D>::t_index m_closest2mate;
        std::vector<t_solve_v2> m_data;
};

template<typename T, unsigned int D>
void Splat<T,D>::strut(void)
{
    t_cells src_cells;
    t_mass src_mass;
    src_cells = m_cells;
    src_mass = m_mass;
    for(unsigned int i = 0; i < m_cells.size(); i++)
    {
        m_cells[i] = m_init;
        m_mass[i] = 0.0;
    }

    Strut<T,D> strt(&src_cells, &src_mass);
    iterate(strt);
}

template<typename T, unsigned int D>
class ReSplat : virtual public Splat<T,D>::Functor
{
    public:
        ReSplat(Splat<T,D> *a_other_splat)
        {
            m_other_splat = a_other_splat;
        }
        virtual ~ReSplat(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            typename Splat<T,D>::t_span span;
            a_splat.spanize(span, a_index);
            m_other_splat->unitSplat(1.0, span, a_splat.get(a_index));
        }
        Splat<T,D> *m_other_splat;
};



template<typename T, unsigned int D>
void Splat<T,D>::resplat(Splat<T,D> *a_other_splat)
{
    ReSplat<T,D> resplat(a_other_splat);
    iterate(resplat);
}

template<typename T, unsigned int D>
class Pull : virtual public Splat<T,D>::Functor
{
    public:
        Pull(Splat<T,D> *a_other_splat)
        {
            m_other_splat = a_other_splat;
        }
        virtual ~Pull(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            typename Splat<T,D>::t_span span;
            a_splat.location(span, a_index);
            a_splat.splat(1.0, span, m_other_splat->sample(span));
        }
        Splat<T,D> *m_other_splat;
};



template<typename T, unsigned int D>
void Splat<T,D>::pull(Splat<T,D> *a_other_splat)
{
    Pull<T,D> p(a_other_splat);
    iterate(p);
}

template<typename T, unsigned int D>
class Push : virtual public Splat<T,D>::Functor
{
    public:
        Push(Splat<T,D> *a_other_splat)
        {
            m_other_splat = a_other_splat;
        }
        virtual ~Push(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            typename Splat<T,D>::t_span span;
            a_splat.location(span, a_index);

            m_other_splat->splat(1.0, span, a_splat.sample(span));
        }
        Splat<T,D> *m_other_splat;
};



template<typename T, unsigned int D>
void Splat<T,D>::push(Splat<T,D> *a_other_splat)
{
    Push<T,D> p(a_other_splat);
    iterate(p);
}

template<typename T, unsigned int D>
class SubGradient : virtual public Splat<T,D>::Functor
{
    public:
        SubGradient(typename Splat<T,D>::t_index &a_index)
        {
            m_index = a_index;
        }
        virtual ~SubGradient(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            unsigned int i_m_cell = a_splat.calc_array_index(m_index);
            typename Splat<T,D>::t_index Dindex = m_index-a_index;
            typename Splat<T,D>::t_span Dspan;
            a_splat.spanize(Dspan, Dindex);
            t_solve_real Dmagnitude = a_splat.magnitude(Dspan);

            if(i_cell == i_m_cell || Dmagnitude < 1.0e-8) { return; }

            typename Splat<T,D>::t_span Dunit;
            for(unsigned int d = 0; d < D; d++)
            {
                Dunit[d] = Dspan[d] / Dmagnitude;
            }




            T value_delta_per_distance = (a_splat.m_cells[i_m_cell] - a_splat.m_cells[i_cell]) / Dmagnitude;

            t_solve_real use_mass = a_splat.m_mass[i_m_cell] / Dmagnitude;
//if(a_splat.m_mass[i_m_cell] > 1.0e-8) { use_mass = 1.0; }
if(a_splat.m_mass[i_m_cell] > 1.0e-8) { use_mass = 1.0 / Dmagnitude; }
else { return; }

            a_splat.gradientSplat(use_mass, a_index, value_delta_per_distance, Dunit);
        }
        typename Splat<T,D>::t_index m_index;
};

template<typename T, unsigned int D>
class Gradient : virtual public Splat<T,D>::Functor
{
    public:
        Gradient()
        {
        }
        virtual ~Gradient(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            SubGradient<T,D> subgradient(a_index);
            a_splat.iterate(subgradient);
        }
};

template<typename T, unsigned int D>
void Splat<T,D>::gradient(void)
{
    Gradient<T,D> gradient_fn;
    iterate(gradient_fn);
}


template<typename T, unsigned int D>
class ReApply : virtual public Splat<T,D>::Functor
{
    public:
        ReApply(typename Splat<T,D>::t_span a_loc, const T &a_value, t_solve_real a_mass)
        {
            m_loc = a_loc;
            m_value = a_value;
            m_mass = a_mass;
        }
        virtual ~ReApply(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            typename Splat<T,D>::t_span loc;
            a_splat.location(loc, a_index);
            t_solve_real sum = 0;
            for(unsigned int i = 0; i < D; i++)
            {
                t_solve_real a = (m_loc[i] - loc[i]) * a_splat.m_scale[i];
                sum += a*a;
            }
            //t_solve_real distance = sum;
            //t_solve_real use_mass = m_mass / (distance + 1.0e-8);
            a_splat.splat(m_mass, loc, m_value);
        }
        typename Splat<T,D>::t_span m_loc;
        T m_value;
        t_solve_real m_mass;
};

template<typename T, unsigned int D>
class Blurt : virtual public Splat<T,D>::Functor
{
    public:
        Blurt(typename Splat<T,D>::t_index &a_index, const T &a_value, t_solve_real a_mass, t_solve_real a_factor)
        {
            m_index = a_index;
            m_value = a_value;
            m_mass = a_mass;
            m_factor = a_factor;
        }
        virtual ~Blurt(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            t_solve_real sum = 0;
            for(unsigned int i = 0; i < D; i++)
            {
                t_solve_real a = (t_solve_real)(m_index[i] - a_index[i]) * a_splat.m_size[i];
                sum += a*a;
            }
            t_solve_real d = fabs(sum);
            t_solve_real metric = m_factor / (m_factor + d);
            if(metric < 0.0) { metric = 0.0; }
            t_solve_real use_mass = 1.0 * m_mass * metric;
            a_splat.massSplat(use_mass, a_index, m_value);
        }
        typename Splat<T,D>::t_index m_index;
        T m_value;
        t_solve_real m_mass;
        t_solve_real m_factor;
};

template<typename T, unsigned int D>
class Blurer : virtual public Splat<T,D>::Functor
{
    public:
        Blurer(typename Splat<T,D>::t_cells *a_cells, typename Splat<T,D>::t_mass *a_mass, t_solve_real a_factor)
        {
            m_src_cells = a_cells;
            m_src_mass = a_mass;
            m_factor = a_factor;
        }
        virtual ~Blurer(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            unsigned int i_cell = a_splat.calc_array_index(a_index);
            T &src_value = (*m_src_cells)[i_cell];
            t_solve_real src_mass = (*m_src_mass)[i_cell];
            {
                Blurt<T,D> blurt(a_index, src_value, src_mass, m_factor);
                a_splat.iterate(blurt);
                //ReApply<T,D> reapply(loc, src_value, src_mass);
                //a_splat.iterate(reapply);
            }
            //fe_fprintf(stderr, "%g\n", a_splat.get(a_index));
            //a_splat.splat(src_mass, loc, src_value);
        }
        typename Splat<T,D>::t_cells *m_src_cells;
        typename Splat<T,D>::t_mass *m_src_mass;
        t_solve_real m_factor;
};

template<typename T, unsigned int D>
void Splat<T,D>::blur(t_solve_real a_factor)
{
    t_cells src_cells;
    t_mass src_mass;
    src_cells = m_cells;
    src_mass = m_mass;
    for(unsigned int i = 0; i < m_cells.size(); i++)
    {
        //m_cells[i] = m_init;
        //m_mass[i] = 0.0;
    }

    Blurer<T,D> blurer(&src_cells, &src_mass, a_factor);
    iterate(blurer);
}

#if 0
template<typename T, unsigned int D>
class Connector : virtual public Splat<T,D>::Functor
{
    public:
        Connector(sp<Scope> spScope, sp<RecordGroup> rg_fill, std::vector<Record> *r_ps, std::vector<Record> *r_cs, std::vector<Record> *r_ks)
        {
            m_spScope = spScope;
            m_rg_fill = rg_fill;

            m_asSolverParticle.bind(spScope);
            m_asParticle.bind(spScope);
            m_asColored.bind(spScope);
            m_asLinearSpring.bind(spScope);

            l_point = m_spScope->declare("l_point");
            m_asSolverParticle.populate(l_point);
            m_asParticle.populate(l_point);
            m_asColored.populate(l_point);

            l_link = m_spScope->declare("l_link");
            m_asLinearSpring.populate(l_link);

            m_pps = r_ps;
            m_pcs = r_cs;
            m_pks = r_ks;
        }
        virtual ~Connector(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            for(unsigned int d = 0; d < D; d++)
            {
                if(a_index[d] < a_splat.m_count[d]-1)
                {
                    unsigned int i_a = a_splat.calc_array_index(a_index);
                    typename Splat<T,D>::t_index other_index;
                    other_index = a_index;
                    other_index[d] += 1;
                    unsigned int i_b = a_splat.calc_array_index(other_index);

                    Record r_link = m_spScope->createRecord(l_link);
                    m_asLinearSpring.left(r_link) = (*m_pps)[i_a];
                    m_asLinearSpring.right(r_link) = (*m_pps)[i_b];
                    m_asLinearSpring.length(r_link) = 0.0;
                    m_asLinearSpring.stiffness(r_link) = 1.0e4;
                    m_asLinearSpring.damping(r_link) = 1.0e1;
                    m_rg_fill->add(r_link);
                }
            }
        }

        std::vector<Record> *m_pps;
        std::vector<Record> *m_pcs;
        std::vector<Record> *m_pks;
        sp<Scope> m_spScope;
        sp<RecordGroup> m_rg_fill;
        AsSolverParticle1D m_asSolverParticle;
        AsParticle1D m_asParticle;
        AsColored m_asColored;
        AsLinearSpringExtended  m_asLinearSpring;
        sp<Layout> l_point;
        sp<Layout> l_link;
};


template<typename T, unsigned int D>
void Splat<T,D>::fillSolve(unsigned int a_aspect)
{
    sp<SemiImplicit1D>  spSemiImplicit;
    sp<Master> spMaster(FLEXI_NEW( gSTMBase.GetMemPool(), gSTMDebug.GetMemTrace(), MWL::Base::BDbgString("Splat") ) Master);
    sp<Registry> spRegistry=spMaster->registry();
    spRegistry->manage("fexCathedralDL");
    sp<Scope> spScope = spRegistry->create("Scope");
    assertMath(spMaster->typeMaster());
    spSemiImplicit = FLEXI_NEW( gSTMBase.GetMemPool(), gSTMDebug.GetMemTrace(), MWL::Base::BDbgString("Splat") ) SemiImplicit1D();
    spSemiImplicit->initialize(spScope);

    {
        sp<SemiImplicit1D::Force> spForce(FLEXI_NEW( gSTMBase.GetMemPool(), gSTMDebug.GetMemTrace(), MWL::Base::BDbgString("Splat") ) LinearSpringForce1D(spSemiImplicit));
        spSemiImplicit->addForce(spForce, true);
    }

    sp<RecordGroup> rg_fill;
    rg_fill = FLEXI_NEW( gSTMBase.GetMemPool(), gSTMDebug.GetMemTrace(), MWL::Base::BDbgString("SetaTireRig") ) RecordGroup();

    AsSolverParticle1D asSolverParticle;
    asSolverParticle.bind(spScope);
    AsParticle1D asParticle;
    asParticle.bind(spScope);
    AsForcePoint1D asForcePoint;
    asForcePoint.bind(spScope);
    AsColored asColored;
    asColored.bind(spScope);

    sp<Layout> l_point = spScope->declare("l_point");
    asSolverParticle.populate(l_point);
    asParticle.populate(l_point);
    asColored.populate(l_point);

    sp<Layout> l_constraint = spScope->declare("l_constraint");
    asForcePoint.populate(l_constraint);
    asColored.populate(l_constraint);

    AsLinearSpringExtended  asLinearSpring;
    asLinearSpring.bind(spScope);
    sp<Layout> l_spring = spScope->declare("l_spring");
    asLinearSpring.populate(l_spring);

    std::vector<Record> r_ps;
    std::vector<Record> r_cs;
    std::vector<Record> r_ks;
    r_ps.resize(m_cells.size());
    r_cs.resize(m_cells.size());
    r_ks.resize(m_cells.size());
    for(unsigned int i = 0; i < m_cells.size(); i++)
    {
        r_ps[i] = spScope->createRecord(l_point);
        asParticle.mass(r_ps[i]) = 1.0;
        asParticle.location(r_ps[i]) = 0.0;
        asParticle.velocity(r_ps[i]) = 0.0;
        rg_fill->add(r_ps[i]);

        r_cs[i] = spScope->createRecord(l_constraint);
        asParticle.location(r_cs[i]) = value(m_cells[i], a_aspect);
        asParticle.velocity(r_cs[i]) = 0.0;
        rg_fill->add(r_cs[i]);

        r_ks[i] = spScope->createRecord(l_spring);
        asLinearSpring.left(r_ks[i]) = r_ps[i];
        asLinearSpring.right(r_ks[i]) = r_cs[i];
        asLinearSpring.length(r_ks[i]) = 0.0;
        asLinearSpring.stiffness(r_ks[i]) = 1.0e8 * m_mass[i];
        asLinearSpring.damping(r_ks[i]) = 1.0e2;
        rg_fill->add(r_ks[i]);
    }

    Connector<T,D> connector(spScope, rg_fill, &r_ps, &r_cs, &r_ks);
    iterate(connector);

    spSemiImplicit->compile(rg_fill);
    std::vector<SemiImplicit1D::Particle> &particles =
        spSemiImplicit->particles();


    t_solve_real dummy;
    for(unsigned int t = 0; t < 1000; t++)
    {
        spSemiImplicit->step(0.001, dummy);
    }

    for(unsigned int i = 0 ; i < r_ps.size(); i++)
    {
        unsigned int i_p = asSolverParticle.index(r_ps[i]);
        fe_fprintf(stderr, "%u %g\n", i_p, particles[i_p].m_location);
        value(m_cells[i], a_aspect) = particles[i_p].m_location;
    }
}
#endif

template<typename T, unsigned int D>
inline typename Splat<T,D>::t_span &Splat<T,D>::safe(
        t_span &a_span) const
{
    for(unsigned int i = 0; i < D; i++)
    {
        if(a_span[i] < 0.0)
        {
            a_span[i] = 0.0;
        }
        else if(a_span[i] > 1.0)
        {
            a_span[i] = 1.0;
        }
    }

    return a_span;
}

template<typename T, unsigned int D>
inline bool Splat<T,D>::safe(t_index &a_index) const
{
    bool ok = true;
    for(unsigned int i = 0; i < D; i++)
    {
        if(!safe(i, a_index)) { ok = false; }
    }
    return ok;
}

template<typename T, unsigned int D>
inline bool Splat<T,D>::safe(unsigned int a_d, t_index &a_index) const
{
    // return true if the requested index is a valid index
    bool rv = true;
    if(m_count[a_d] == 1)
    {
        if(a_index[a_d] != 0)
        {
            rv = false;
            a_index[a_d] = 0;
        }
    }
    else if(m_count[a_d] == 2)
    {
        if(a_index[a_d] < 0)
        {
            rv = false;
            a_index[a_d] = 0;
        }
        else if(a_index[a_d] > 1)
        {
            rv = false;
            a_index[a_d] = 1;
        }
    }
    else
    {
        if(a_index[a_d] < 0)
        {
            rv = false;
            a_index[a_d] = 0;
        }
        else if(a_index[a_d] >= m_count[a_d])
        {
            rv = false;
            a_index[a_d] = m_count[a_d]-1;
        }
    }

    return rv;
}

template<typename T, unsigned int D>
inline void Splat<T,D>::create(const t_index &a_count, const T &a_init)
{
    m_totalCount = 1;
    for(unsigned int i = 0; i < D; i++)
    {
        m_count[i] = a_count[i];
        m_size[i] = 1.0/(m_count[i]-1);
        m_totalCount *= m_count[i];
    }

    m_init = a_init;
    m_cells.resize(m_totalCount, a_init);
    m_mass.resize(m_totalCount, 0.0);
    m_gradients.resize(m_totalCount);
    m_gradient_mass.resize(m_totalCount, 0.0);
    for(unsigned int i = 0; i < m_totalCount; i++)
    {
        for(unsigned int d = 0; d < D; d++)
        {
            zero(m_gradients[i].m_gradient[d]);
        }
    }
}

template<typename T, unsigned int D>
void Splat<T,D>::clone(Splat<T,D> *a_other_splat)
{
    m_totalCount = a_other_splat->m_totalCount;
    for(unsigned int i = 0; i < D; i++)
    {
        m_count[i] = a_other_splat->m_count[i];
        m_size[i] = a_other_splat->m_size[i];
    }
    m_init = a_other_splat->m_init;
    m_cells.resize(m_totalCount, m_init);
    //m_mass.resize(m_totalCount, 0.0);
    //m_gradients.resize(m_totalCount);
    //m_gradient_mass.resize(m_totalCount, 0.0);
    for(unsigned int i = 0; i < m_totalCount; i++)
    {
        for(unsigned int d = 0; d < D; d++)
        {
            //zero(m_gradients[i].m_gradient[d]);
            m_cells[i][d] = a_other_splat->m_cells[i][d];
        }
    }
    for(unsigned int d = 0; d < D; d++)
    {
        m_root[d] = a_other_splat->m_root[d];
        m_scale[d] = a_other_splat->m_scale[d];
    }
}


template<typename T, unsigned int D>
inline const T &Splat<T,D>::get(const t_index &a_index) const
{
    t_index safe_index(a_index);
    safe(safe_index);
    return m_cells[calc_array_index(safe_index)];
}

template<typename T, unsigned int D>
inline const T &Splat<T,D>::direct(const t_index &a_index) const
{
    //t_index safe_index(a_index);
    //safe(safe_index);
    //return m_cells[calc_array_index(safe_index)];
    return m_cells[calc_array_index(a_index)];
}

template<typename T, unsigned int D>
inline const t_solve_real &Splat<T,D>::getMass(const t_index &a_index) const
{
    t_index safe_index(a_index);
    safe(safe_index);
    return m_mass[calc_array_index(safe_index)];
}

template<typename T, unsigned int D>
inline void Splat<T,D>::set(const t_index &a_index, const T &a_value)
{
    t_index safe_index(a_index);
    safe(safe_index);
    m_cells[calc_array_index(safe_index)] = a_value;
}

template<typename T, unsigned int D>
inline void Splat<T,D>::setMass(const t_index &a_index, const t_solve_real &a_value)
{
    t_index safe_index(a_index);
    safe(safe_index);
    m_mass[calc_array_index(safe_index)] = a_value;
}

template<typename T, unsigned int D>
inline unsigned int Splat<T,D>::calc_array_index(
        const t_index &a_index) const
{
    unsigned int index = 0;
    unsigned int step = 1;
    for(int d = D-1; d >= 0; d--)
    {
        index += a_index[d] * step;
        step *= m_count[d];
    }
    return index;
}

template<typename T, unsigned int D>
inline void Splat<T,D>::iterate(Functor &a_functor)
{
    t_index index;
    rIterate(D-1, index, a_functor);
}

template<typename T, unsigned int D>
inline void Splat<T,D>::rIterate(unsigned int a_d, t_index &a_index,
        Functor &a_functor)
{
    for(int i = 0; i < count()[a_d]; i++)
    {
        a_index[a_d] = i;
        if(a_d == 0)
        {
            a_functor.process(*this, a_index);
        }
        else
        {
            rIterate(a_d-1, a_index, a_functor);
        }
    }
}

template<typename T, unsigned int D>
class Dumper : virtual public Splat<T,D>::Functor
{
    public:
        Dumper(const char *a_filename)
        {
            m_fp = fopen(a_filename, "w");
        }
        virtual ~Dumper(void)
        {
            fclose(m_fp);
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            if(D <=2 || a_index[2] == 2)
            {
            //fe_fprintf(m_fp, "%g %g %g\n", (Real)(a_index[0]), (Real)(a_index[1]), a_splat.get(a_index)[1]);
                typename Splat<T,D>::t_span loc;
                a_splat.location(loc, a_index);
                fe_fprintf(m_fp, "%g %g %g\n", loc[0], loc[1], a_splat.get(a_index)[0]);
            //fe_fprintf(m_fp, "%g %g %g\n", (Real)(a_index[0]), (Real)(a_index[1]), a_splat.getMass(a_index));
            }
        }

        FILE *m_fp;
};


template<typename T, unsigned int D>
inline void Splat<T,D>::dump(const char *a_filename)
{
    Dumper<T,D> dumper(a_filename);
    iterate(dumper);
}

typedef Splat<t_solve_v3,2>     t_image_splat;


template<typename T, unsigned int D>
class DBSaver : virtual public Splat<T,D>::Functor
{
    public:
        DBSaver(const char *a_filename)
        {
            m_fp = fopen(a_filename, "w");
        }
        virtual ~DBSaver(void)
        {
            fclose(m_fp);
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            fe_fprintf(m_fp, "%g %g ", a_splat.get(a_index)[0], a_splat.get(a_index)[1]);
        }

        FILE *m_fp;
};


template<typename T, unsigned int D>
inline void Splat<T,D>::dbsave(const char *a_filename)
{
    DBSaver<T,D> saver(a_filename);
    iterate(saver);
}

#if 0
template<typename T, unsigned int D>
class Saver : virtual public Splat<T,D>::Functor
{
    public:
        Saver(Record &r_splat)
        {
            m_r_splat = r_splat;
        }
        virtual ~Saver(void)
        {
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
        }

        Record m_r_splat;
};
#endif

#if 0
inline void save(Splat<t_solve_v2,2> &a_splat, Record &r_splat)
{
    AsSplat22 asSplat;
    sp<Scope> spScope = r_splat.layout()->scope();

    asSplat.bind(spScope);

    asSplat.count(r_splat) = Vector2i(a_splat.m_count[0], a_splat.m_count[1]);
}
#endif

template<typename T, unsigned int D>
class Loader : virtual public Splat<T,D>::Functor
{
    public:
        Loader(const char *a_filename)
        {
            m_fp = fopen(a_filename, "r");
        }
        virtual ~Loader(void)
        {
            fclose(m_fp);
        }
        virtual void process(Splat<T,D> &a_splat, typename Splat<T,D>::t_index &a_index)
        {
            float f0,f1;
            fscanf(m_fp, "%g %g", &f0, &f1);
            a_splat.set(a_index, t_solve_v2(f0,f1));
        }

        FILE *m_fp;
};


#if 0
Requires new state access methodology...back to tires first

template<typename T, unsigned int D>
inline void Splat<T,D>::dbload(const char *a_filename)
{
    Loader<T,D> loader(a_filename);
    iterate(loader);
}

typedef Splat<t_solve_v3,2>     t_image_splat;

class TestCase
{
    public:
        TestCase(sp<SetaTire> a_tire, const String &a_label, const String &a_testname, const String &a_id);
        TestCase(sp<TireManagerI> a_tm, const String &a_label, const String &a_testname, const String &a_id);
virtual ~TestCase(void);
virtual void run(void);
virtual void runCase(void) = 0;

        void mark(t_tire_real a_x, t_tire_real a_y);
        void line(t_tire_real a_x, t_tire_real a_y);
        void newline(void);
        void annotate(t_tire_real a_x, t_tire_real a_y, const String &a_text);
        void annotate(t_tire_real a_x, t_tire_real a_y, t_tire_real a_value);
        void hilight(t_scalar_tire_data a_datum, t_tire_real a_center, t_tire_real a_falloff, t_tire_v3 a_color);
        void hline(t_tire_real a_x, t_tire_real a_y, const t_tire_v3 &a_color, const String &a_text);
        void hline(t_tire_real a_x, t_tire_real a_y, const t_tire_v3 &a_color, t_tire_real a_value);
        void hline(t_tire_real a_x, t_tire_real a_y, const t_tire_v3 &a_color);
        void vline(t_tire_real a_x, t_tire_real a_y, const t_tire_v3 &a_color, const String &a_text);
        void vline(t_tire_real a_x, t_tire_real a_y, const t_tire_v3 &a_color, t_tire_real a_value);
        void vline(t_tire_real a_x, t_tire_real a_y, const t_tire_v3 &a_color);
        void bar(t_tire_real a_x, t_tire_real a_y, const t_tire_v3 &a_color);
        void image(const char *a_filename);
        void dumpVis(void);

        const String &testname(void) { return m_testname; }
        const String &id(void) { return m_id; }

        class Hilight
        {
            public:
                Hilight(t_scalar_tire_data a_datum, t_tire_real a_center, t_tire_real a_falloff, t_tire_v3 a_color)
                {
                    m_datum = a_datum;
                    m_color = a_color;
                    m_center = a_center;
                    m_falloff = a_falloff;
                };
                t_scalar_tire_data m_datum;
                t_tire_v3 m_color;
                t_tire_real m_center;
                t_tire_real m_falloff;
        };

        class Annotation
        {
            public:
                Annotation(t_tire_real a_x, t_tire_real a_y, const String &a_text) :
                    m_x(a_x), m_y(a_y), m_text(a_text) {}
                t_tire_real m_x;
                t_tire_real m_y;
                String      m_text;
        };
    protected:
        sp<SetaTire>                m_spTire;
        t_image_splat               m_splat;
        SetaTireRig                 m_rig;
        t_stdvector< Hilight >      m_hilights;
        String                      m_label;
        String                      m_testname;
        String                      m_id;
        t_stdvector< Annotation >   m_annotations;
        bool                        m_line_started;
        t_image_splat::t_span       m_loc;
        t_tire_v3                   m_default_color;
        String                      m_data_path;
        unsigned int                m_frame;
};
#endif


} /* namespace ext */
} /* namespace fe */


#endif /* __solve_Splat_h__ */