#ifndef SOFA_BRANCHINGIMAGE_ImageValuesFromPositions_H
#define SOFA_BRANCHINGIMAGE_ImageValuesFromPositions_H

#include <BranchingImage/config.h>
#include <image/ImageValuesFromPositions.h>
#include "BranchingImage.h"

namespace sofa
{
namespace component
{
namespace engine
{


/**
 * Get image intensities at sample locations
 */

/// Specialization for branching Image
template <class T>
struct ImageValuesFromPositionsSpecialization<defaulttype::BranchingImage<T>>
{
    typedef ImageValuesFromPositions<defaulttype::BranchingImage<T>> ImageValuesFromPositionsT;

    static void update(ImageValuesFromPositionsT& This)
    {
        typedef typename ImageValuesFromPositionsT::Real Real;
        typedef typename ImageValuesFromPositionsT::Coord Coord;
        typedef typename ImageValuesFromPositionsT::ImageTypes::VoxelIndex VoxelIndex;

        typename ImageValuesFromPositionsT::raTransform inT(This.transform);
        typename ImageValuesFromPositionsT::raPositions pos(This.position);

        typename ImageValuesFromPositionsT::raImage in(This.image);
        if(in->isEmpty()) return;
        const typename ImageValuesFromPositionsT::ImageTypes& img = in.ref();

        typename ImageValuesFromPositionsT::waValues val(This.values);
        Real outval=This.outValue.getValue();
        val.resize(pos.size());


        switch(This.Interpolation.getValue().getSelectedId())
        {
        case INTERPOLATION_CUBIC :
            This.serr<<"Cubic Interpolation not implemented for branching images -> nearest interpolation"<<This.sendl;
//        {
//            for(unsigned int i=0; i<pos.size(); i++)
//            {
//                Coord Tp = inT->toImage(pos[i]);
//                if(!in->isInside(Tp[0],Tp[1],Tp[2]))  val[i] = outval;
//                else val[i] = (Real)img.cubic_atXYZ(Tp[0],Tp[1],Tp[2],0,(T)outval,cimg_library::cimg::type<T>::min(),cimg_library::cimg::type<T>::max());
//            }
//        }
//            break;

        case INTERPOLATION_LINEAR :
            This.serr<<"Linear Interpolation not implemented for branching images -> nearest interpolation"<<This.sendl;
//        {
//            for(unsigned int i=0; i<pos.size(); i++)
//            {
//                Coord Tp = inT->toImage(pos[i]);
//                if(!in->isInside(Tp[0],Tp[1],Tp[2])) val[i] = outval;
//                else val[i] = (Real)img.linear_atXYZ(Tp[0],Tp[1],Tp[2],0,(T)outval);
//            }
//        }
//            break;

        default : // NEAREST
        {
            for(unsigned int i=0; i<pos.size(); i++)
            {
                Coord Tp = inT->toImageInt(pos[i]);
                if(!in->isInside((int)Tp[0],(int)Tp[1],(int)Tp[2]))  val[i] = outval;
                else
                {
                    const VoxelIndex vi (in->index3Dto1D(Tp[0],Tp[1],Tp[2]), 0);  // take first superimposed voxel
                    if(in->Nb_superimposed(vi.index1d, This.time)) val[i] = (Real)img(vi,0,This.time);  // take first channel
                    else val[i] = outval;
                }
            }
        }
            break;
        }

    }

};


} // namespace engine
} // namespace component
} // namespace sofa

#endif // SOFA_BRANCHINGIMAGE_ImageValuesFromPositions_H
