Subversion Repositories gelsvn

/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */

#include "GradientFilter.h"
#include "TrilinFilter.h"
#include "Neighbours.h"


using namespace CGLA;

namespace Geometry
{
        template<class GridT>
        bool TrilinFilter<GridT>::in_domain(const CGLA::Vec3f& v) const 
        {
                Vec3i c0i(v);
                return grid->in_domain(c0i) && grid->in_domain(c0i+Vec3i(1));
        }

        template<class GridT>
        float TrilinFilter<GridT>::operator()(const CGLA::Vec3f& v) const
        {
                Vec3i c0i(v);

                const float alpha = v[0] - float(c0i[0]);
                const float beta  = v[1] - float(c0i[1]);
                const float gamm  = v[2] - float(c0i[2]);

                float m_alpha = 1.0 - alpha;
                float m_beta  = 1.0 - beta;
                float m_gamm  = 1.0 - gamm;

                float voxels[8];
                for(int i=0;i<8;++i)
                        voxels[i] = (*grid)[c0i+Geometry::CubeCorners8i[i]];
  
                float f = 0.0f;
                f += (m_alpha*m_beta*m_gamm)*voxels[0];
                f += (alpha*m_beta*m_gamm)*voxels[1];
                f += (m_alpha*beta*m_gamm)*voxels[2];
                f += (alpha*beta*m_gamm)*voxels[3];
                f += (m_alpha*m_beta*gamm)*voxels[4];
                f += (alpha*m_beta*gamm)*voxels[5];
                f += (m_alpha*beta*gamm)*voxels[6];
                f += (alpha*beta*gamm)*voxels[7];
                
                return f;
        }

        template<class GridT>
        Vec3f TrilinFilter<GridT>::grad(const CGLA::Vec3f& v) const
        {
                Vec3i c0i(v);

                const float alpha = v[0] - float(c0i[0]);
                const float beta  = v[1] - float(c0i[1]);
                const float gamm  = v[2] - float(c0i[2]);

                float m_alpha = 1.0 - alpha;
                float m_beta  = 1.0 - beta;
                float m_gamm  = 1.0 - gamm;

                float voxels[8];
                for(int i=0;i<8;++i)
                        voxels[i] = (*grid)[c0i+Geometry::CubeCorners8i[i]];
  
                float fx = 0.0f;
                fx += -m_beta*m_gamm*float(voxels[0]);
                fx += m_beta*m_gamm*float(voxels[1]);
                fx += -beta*m_gamm*float(voxels[2]);
                fx += beta*m_gamm*float(voxels[3]);
                fx += -m_beta*gamm*float(voxels[4]);
                fx += m_beta*gamm*float(voxels[5]);
                fx += -beta*gamm*float(voxels[6]);
                fx += beta*gamm*float(voxels[7]);

                float fy = 0.0f;
                fy += -m_alpha*m_gamm*float(voxels[0]);
                fy += -alpha*m_gamm*float(voxels[1]);
                fy += m_alpha*m_gamm*float(voxels[2]);
                fy += alpha*m_gamm*float(voxels[3]);
                fy += -m_alpha*gamm*float(voxels[4]);
                fy += -alpha*gamm*float(voxels[5]);
                fy += m_alpha*gamm*float(voxels[6]);
                fy += alpha*gamm*float(voxels[7]);

                float fz = 0.0f;
                fz += -m_alpha*m_beta*float(voxels[0]);
                fz += -alpha*m_beta*float(voxels[1]);
                fz += -m_alpha*beta*float(voxels[2]);
                fz += -alpha*beta*float(voxels[3]);
                fz += m_alpha*m_beta*float(voxels[4]);
                fz += alpha*m_beta*float(voxels[5]);
                fz += m_alpha*beta*float(voxels[6]);
                fz += alpha*beta*float(voxels[7]);

                return Vec3f(fx, fy, fz);
        }

                template class  TrilinFilter<RGrid<float> >;
                template class  TrilinFilter<RGrid<unsigned char> >;
                template class  TrilinFilter<RGrid<unsigned short> >;
                template class  TrilinFilter<RGrid<short> >;
}