#ifndef SOFA_BRANCHINGIMAGE_BranchingCellOffsetsFromConnectLabels_H
#define SOFA_BRANCHINGIMAGE_BranchingCellOffsetsFromConnectLabels_H

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

#include <sofa/core/objectmodel/Event.h>
#include <sofa/simulation/AnimateEndEvent.h>

#include <sofa/defaulttype/Vec.h>
#include <sofa/helper/SVector.h>

namespace sofa
{
namespace component
{
namespace engine
{

/**
 * Returns offsets of superimposed voxels at positions where point labels and image labels overlap
 */


template <class _ImageTypes>
class BranchingCellOffsetsFromConnectLabels : public core::DataEngine
{
public:
    typedef core::DataEngine Inherited;
    SOFA_CLASS(SOFA_TEMPLATE(BranchingCellOffsetsFromConnectLabels,_ImageTypes),Inherited);

    typedef SReal Real;

    //@name Image data
    /**@{*/
    typedef _ImageTypes ImageTypes;
    typedef typename ImageTypes::T T;
    typedef helper::ReadAccessor<Data< ImageTypes > > raImage;
    Data< ImageTypes > d_image;
    /**@}*/

    //@name Transform data
    /**@{*/
    typedef defaulttype::ImageLPTransform<Real> TransformType;
    typedef typename TransformType::Coord Coord;
    typedef helper::ReadAccessor<Data< TransformType > > raTransform;
    Data< TransformType > d_transform;
    /**@}*/

    //@name Label data
    /**@{*/
    typedef T LabelType;
    typedef helper::vector<T> Vlabels;
    typedef helper::vector< helper::SVector<T> > VecVlabels;
    typedef helper::ReadAccessor<Data< Vlabels > > raVLabels;
    typedef helper::ReadAccessor<Data< VecVlabels > > raVecVLabels;
    Data< VecVlabels > d_pointLabels;
    Data< Vlabels > d_connectLabels;
    /**@}*/

    //@name position data
    /**@{*/
    typedef helper::vector<defaulttype::Vec<3,Real> > SeqPositions;
    typedef helper::ReadAccessor<Data< SeqPositions > > raPositions;
    Data< SeqPositions > d_position;
    /**@}*/

    //@name parameters
    /**@{*/
    Data<bool> d_useGlobalIndices;
    Data<bool> d_trim;
    /**@}*/

    //@name outputs
    /**@{*/
    typedef int OffsetType;
    typedef helper::vector<OffsetType> VOffset;
    typedef helper::WriteOnlyAccessor<Data< VOffset > > waOffset;
    Data< VOffset > d_cell;
    typedef unsigned int IndexType;
    typedef helper::vector<IndexType> VIndex;
    typedef helper::WriteOnlyAccessor<Data< VIndex > > waIndex;
    Data< VIndex > d_indices;
    /**@}*/

    virtual std::string getTemplateName() const    { return templateName(this);    }
    static std::string templateName(const BranchingCellOffsetsFromConnectLabels<ImageTypes>* = NULL) { return ImageTypes::Name();    }

    BranchingCellOffsetsFromConnectLabels()    :   Inherited()
      , d_image(initData(&d_image,ImageTypes(),"image",""))
      , d_transform(initData(&d_transform,TransformType(),"transform",""))
      , d_pointLabels(initData(&d_pointLabels,"pointLabels","lists of labels associated to each point"))
      , d_connectLabels(initData(&d_connectLabels,"connectLabels","Pairs of label to be connected accross different input images"))
      , d_position(initData(&d_position,SeqPositions(),"position","input positions"))
      , d_useGlobalIndices(initData(&d_useGlobalIndices,false,"useGlobalIndices","if true, global cell indices are computed"))
      , d_trim(initData(&d_trim,true,"trim","if true, returns ouput only for connected points"))
      , d_cell( initData ( &d_cell,"cell","cell offsets" ) )
      , d_indices( initData ( &d_indices,"indices","indices of connected points" ) )
      , time((unsigned int)0)
    {
        this->addAlias(&d_transform, "branchingTransform");
        this->addAlias(&d_image, "branchingImage");
        d_transform.setReadOnly(true);
        d_image.setReadOnly(true);

        f_listening.setValue(true);
    }

    virtual void init()
    {
        addInput(&d_image);
        addInput(&d_transform);
        addInput(&d_pointLabels);
        addInput(&d_connectLabels);
        addInput(&d_position);
        addInput(&d_useGlobalIndices);
        addInput(&d_trim);
        addOutput(&d_cell);
        addOutput(&d_indices);
        setDirtyValue();
    }

    virtual void reinit() { update(); }

protected:

    unsigned int time;

    virtual void update()
    {
        updateAllInputsIfDirty(); // easy to ensure that all inputs are up-to-date

        cleanDirty();

        raImage in(d_image);
        raTransform inT(d_transform);
        raPositions pos(d_position);
        raVLabels connectLabels(d_connectLabels);
        raVecVLabels pointLabels(d_pointLabels);
        bool glob=d_useGlobalIndices.getValue();
        bool trim=d_trim.getValue();

        if(pos.size()!=pointLabels.size())
        {
            serr<<"position and pointLabels data have different sizes"<<sendl;
            return;
        }

        // precompute set for fast look up
        typedef std::pair<T,T>  TPair;
        std::set<TPair> connectS;
        for(unsigned int i=0;i<connectLabels.size()/2;i++)
        {
            connectS.insert(TPair(connectLabels[2*i],connectLabels[2*i+1]));
            connectS.insert(TPair(connectLabels[2*i+1],connectLabels[2*i]));
        }

        // get images at time t
        const typename ImageTypes::BranchingImage3D& img = in->imgList[this->time];

        waOffset cell(d_cell); cell.clear();
        waIndex indices(d_indices); indices.clear();

        OffsetType outOffset=-1;

        for(size_t i=0; i<pos.size(); i++)
        {
            OffsetType offset=outOffset;
            Coord Tp = inT->toImageInt(pos[i]);
            if(in->isInside((int)Tp[0],(int)Tp[1],(int)Tp[2]))
            {
                typename ImageTypes::VoxelIndex vi (in->index3Dto1D(Tp[0],Tp[1],Tp[2]), 0);
                bool found=false;
                for(vi.offset = 0 ; vi.offset<img[vi.index1d].size() ; vi.offset++)
                {
                    T v= img[vi.index1d][vi.offset][0]; // assume that labels are stored in first channel
                    for(size_t j=0;j<pointLabels[i].size();++j)
                        if(connectS.find(TPair(v,pointLabels[i][j]))!=connectS.end())
                        {
                            found=true;
                            break;
                        }
                    if(found) break;
                }
                if(found)
                {
                    offset=vi.offset+1;
                    if(glob) {for(unsigned int i1d=0;i1d<vi.index1d;i1d++) offset+=img[i1d].size(); offset-=1;}
                }
            }
            if(!trim || offset!=outOffset)
            {
                cell.push_back(offset);
                indices.push_back(i);
            }
        }
    }


    void handleEvent(sofa::core::objectmodel::Event *event)
    {
        if (simulation::AnimateEndEvent::checkEventType(event))
        {
            raImage in(this->d_image);
            raTransform inT(this->d_transform);

            // get current time modulo dimt
            const unsigned int dimt=in->getDimensions()[4];
            if(!dimt) return;
            Real t=inT->toImage(this->getContext()->getTime()) ;
            t-=(Real)((int)((int)t/dimt)*dimt);
            t=(t-floor(t)>0.5)?ceil(t):floor(t); // nearest
            if(t<0) t=0.0; else if(t>=(Real)dimt) t=(Real)dimt-1.0; // clamp

            if(this->time!=(unsigned int)t) { this->time=(unsigned int)t; update(); }
        }
    }

};


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

#endif // SOFA_BRANCHINGIMAGE_BranchingCellOffsetsFromConnectLabels_H
