DiscreteExteriorCalculus.ih

/**
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 **/
 
/**
 * @file DiscreteExteriorCalculus.ih
 * @author Pierre Gueth (\c pierre.gueth@liris.cnrs.fr )
 * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
 *
 * @date 2014/03/27
 *
 * Implementation of inline methods defined in DiscreteExteriorCalculus.h
 *
 * This file is part of the DGtal library.
 */
 
///////////////////////////////////////////////////////////////////////////////
// IMPLEMENTATION of inline methods.
///////////////////////////////////////////////////////////////////////////////
 
template <DGtal::Dimension dim, typename TInteger>
size_t
DGtal::hash_value(const DGtal::KhalimskyCell<dim, TInteger>& cell)
{
    return std::hash<size_t>()( boost::hash_range( cell.preCell().coordinates.begin(), cell.preCell().coordinates.end()) );
}
 
///////////////////////////////////////////////////////////////////////////////
// ----------------------- Standard services ------------------------------
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::DiscreteExteriorCalculus()
    : myKSpace(), myCachedOperatorsNeedUpdate(true), myIndexesNeedUpdate(false)
{
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <typename TDomain>
void
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::initKSpace(DGtal::ConstAlias<TDomain> _domain)
{
    BOOST_CONCEPT_ASSERT(( concepts::CDomain<TDomain> ));
 
    // FIXME borders are removed from set => better not initialize kspace
    // FIXME should be open or closed? => closed = true
    const bool kspace_init_ok = const_cast<KSpace&>(myKSpace).init(_domain->lowerBound(), _domain->upperBound(), true);
    ASSERT(kspace_init_ok);
    boost::ignore_unused_variable_warning(kspace_init_ok);
}
 
///////////////////////////////////////////////////////////////////////////////
// Interface - public :
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
bool
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::eraseCell(const Cell& _cell)
{
    typename Properties::iterator iter_property = myCellProperties.find(_cell);
    if (iter_property == myCellProperties.end())
        return false;
 
    myCellProperties.erase(iter_property);
 
    myIndexesNeedUpdate = true;
    myCachedOperatorsNeedUpdate = true;
 
    return true;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
bool
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::insertSCell(const SCell& signed_cell)
{
    return insertSCell(signed_cell, 1, 1);
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
bool
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::insertSCell(const SCell& signed_cell, const Scalar& primal_size, const Scalar& dual_size)
{
    const Cell cell = myKSpace.unsigns(signed_cell);
    const DGtal::Dimension cell_dim = myKSpace.uDim(cell);
    Property property;
    property.primal_size = primal_size;
    property.dual_size = dual_size;
    property.index = std::numeric_limits<Index>::max();
    property.flipped = ( myKSpace.sSign(signed_cell) == KSpace::NEG );
 
    if (cell_dim == 0) ASSERT_MSG( property.primal_size == 1, "0-cell primal size must be equal to 1" );
    if (cell_dim == dimEmbedded) ASSERT_MSG( property.dual_size == 1, "n-cell dual size must be equal to 1" );
    ASSERT_MSG( cell_dim <= dimEmbedded, "wrong cell dimension" );
    ASSERT_MSG( cell_dim != 0 || !property.flipped , "can't insert negative 0-cells" );
    ASSERT_MSG( cell_dim != dimAmbient || !property.flipped , "can't insert negative n-cells" );
 
    std::pair<typename Properties::iterator, bool> insert_pair = myCellProperties.insert(std::make_pair(cell, property));
    if (!insert_pair.second) insert_pair.first->second = property;
 
    ASSERT( insert_pair.first->first == cell );
    ASSERT( insert_pair.first->second.dual_size == property.dual_size );
    ASSERT( insert_pair.first->second.primal_size == property.primal_size );
    ASSERT( insert_pair.first->second.index == property.index );
    ASSERT( insert_pair.first->second.flipped == property.flipped );
 
    myIndexesNeedUpdate = true;
    myCachedOperatorsNeedUpdate = true;
 
    return insert_pair.second;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
void
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::resetSizes()
{
    for (typename Properties::iterator pi=myCellProperties.begin(), pe=myCellProperties.end(); pi!=pe; pi++)
    {
        pi->second.primal_size = 1;
        pi->second.dual_size = 1;
    }
 
    myCachedOperatorsNeedUpdate = true;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
std::string
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::className() const
{
    return "Calculus";
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
void
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::selfDisplay(std::ostream & os) const
{
    os << "[dec";
    for (DGtal::Order order=0; order<=dimEmbedded; order++)
        os << " | primal " << order << "-cells <-> dual " << dimEmbedded-order << "-cells (" << kFormLength(order, PRIMAL) << ")";
    os << "]";
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Order order, DGtal::Duality duality, typename TConstIterator>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::reorder(const TConstIterator& begin_range, const TConstIterator& end_range) const
{
    BOOST_STATIC_ASSERT(( boost::is_convertible<typename TConstIterator::value_type, const SCell>::value ));
 
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
 
    typedef typename TLinearAlgebraBackend::Triplet Triplet;
    typedef std::vector<Triplet> Triplets;
    Triplets triplets;
 
    Index original_index = 0;
    for (TConstIterator ci=begin_range; ci!=end_range; ci++)
    {
        const SCell signed_cell = *ci;
        const Dimension dim = myKSpace.sDim(signed_cell);
        if (dim != actualOrder(order, duality)) continue;
        const Cell cell = myKSpace.unsigns(signed_cell);
        const Index calculus_index = getCellIndex(cell);
        triplets.push_back( Triplet(original_index, calculus_index, 1) );
        original_index++;
    }
 
    const Index length = kFormLength(order, duality);
    ASSERT( triplets.size() == static_cast<typename Triplets::size_type>(length) );
    SparseMatrix reorder_matrix(length, length);
    reorder_matrix.setFromTriplets(triplets.begin(), triplets.end());
 
    typedef LinearOperator<Self, order, duality, order, duality> ReorderOperator;
    return ReorderOperator(*this, reorder_matrix);
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Order order, DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::identity() const
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
 
    typedef LinearOperator<Self, order, duality, order, duality> Operator;
    Operator id(*this);
    id.myContainer.setIdentity();
    return id;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 0, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::laplace() const
{
    typedef DGtal::LinearOperator<Self, 0, duality, 1, duality> Derivative;
    typedef DGtal::LinearOperator<Self, 1, duality, 0, duality> Antiderivative;
    const Derivative d = derivative<0, duality>();
    const Antiderivative ad = antiderivative<1, duality>();
    return ad * d;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 0, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::heatLaplace(const typename DenseVector::Scalar& h, 
  const typename DenseVector::Scalar& t, const typename DenseVector::Scalar& K) const
{
  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
  
  typedef typename TLinearAlgebraBackend::Triplet Triplet;
  typedef LinearOperator<Self, 0, duality, 0, duality> Operator;
 
  const typename DenseVector::Scalar cut = K * sqrt(2. * t);
 
  DGtal::CanonicSCellEmbedder<KSpace> canonicSCellEmbedder(myKSpace);
 
  typedef std::vector<Triplet> Triplets;
  Triplets triplets;
 
  for (Index i = 0; i < kFormLength(0, duality); i++)
  {
    trace.progressBar( i, kFormLength( 0, duality ) - 1 );
 
    const SCell signed_cell_i = myIndexSignedCells[ actualOrder(0, duality) ][i];
    const auto p_i = double( h ) * canonicSCellEmbedder( signed_cell_i );
 
    typename DenseVector::Scalar total_weight = 0.;
 
    for(Index j = 0; j < kFormLength(0, duality); j++)
    {
      if(i == j) continue;
 
      const SCell signed_cell_j = myIndexSignedCells[ actualOrder(0, duality) ][j];
      const auto p_j = double( h ) * canonicSCellEmbedder( signed_cell_j );
 
      const typename DenseVector::Scalar l2_distance = (p_i - p_j).norm();
      if(l2_distance < cut)
      {
        const typename Properties::const_iterator iter_property = myCellProperties.find( myKSpace.unsigns( signed_cell_j ) );
        if ( iter_property == myCellProperties.end() )
          continue;
 
        const typename DenseVector::Scalar measure = (duality == DUAL) ? iter_property->second.primal_size : iter_property->second.dual_size;
        const typename DenseVector::Scalar laplace_value = measure * exp(- l2_distance * l2_distance / (4. * t)) * ( 1. / (t * pow(4. * M_PI * t, dimEmbedded / 2.)) );
 
        triplets.push_back( Triplet(i, j, laplace_value) );
        total_weight -= laplace_value;
      }
    }
 
    //_operator.myContainer.insert( i, i ) = total_weight;
    triplets.push_back( Triplet( i, i, total_weight ) );
  }
  
  Operator _operator( *this );
  ASSERT( _operator.myContainer.rows() == kFormLength(0, duality) );
  ASSERT( _operator.myContainer.cols() == kFormLength(0, duality) );
  _operator.myContainer.setFromTriplets(triplets.begin(), triplets.end());
  _operator.myContainer.makeCompressed();
 
  return _operator;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Order order, DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order-1, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::antiderivative() const
{
    BOOST_STATIC_ASSERT(( order > 0 ));
    BOOST_STATIC_ASSERT(( order <= dimEmbedded ));
 
    typedef DGtal::LinearOperator<Self, order, duality, dimEmbedded-order, OppositeDuality<duality>::duality> FirstHodge;
    typedef DGtal::LinearOperator<Self, dimEmbedded-order, OppositeDuality<duality>::duality, dimEmbedded-order+1, OppositeDuality<duality>::duality> Derivative;
    typedef DGtal::LinearOperator<Self, dimEmbedded-order+1, OppositeDuality<duality>::duality, order-1, duality> SecondHodge;
    const FirstHodge h_first = hodge<order, duality>();
    const Derivative d = derivative<dimEmbedded-order, OppositeDuality<duality>::duality>();
    const SecondHodge h_second = hodge<dimEmbedded-order+1, OppositeDuality<duality>::duality>();
    const Scalar sign = ( order*(dimEmbedded-order)%2 == 0 ? 1 : -1 );
    return sign * h_second * d * h_first;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Order order, DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order+1, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::derivative() const
{
    BOOST_STATIC_ASSERT(( order >= 0 ));
    BOOST_STATIC_ASSERT(( order < dimEmbedded ));
 
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
 
    typedef typename TLinearAlgebraBackend::Triplet Triplet;
    typedef std::vector<Triplet> Triplets;
    Triplets triplets;
 
    // iterate over output form values
    for (Index index_output=0; index_output<kFormLength(order+1, duality); index_output++)
    {
        const SCell signed_cell = myIndexSignedCells[actualOrder(order+1, duality)][index_output];
 
        // find cell border
        typedef typename KSpace::SCells Border;
        const Border border = ( duality == PRIMAL ? myKSpace.sLowerIncident(signed_cell) : myKSpace.sUpperIncident(signed_cell) );
 
        // iterate over cell border
        for (typename Border::const_iterator bi=border.begin(), bie=border.end(); bi!=bie; bi++)
        {
            const SCell signed_cell_border = *bi;
            ASSERT( myKSpace.sDim(signed_cell_border) == actualOrder(order, duality) );
 
            const typename Properties::const_iterator iter_property = myCellProperties.find(myKSpace.unsigns(signed_cell_border));
            if ( iter_property == myCellProperties.end() )
                continue;
 
            const Index index_input = iter_property->second.index;
            ASSERT( index_input < kFormLength(order, duality) );
 
            const bool flipped_border = ( myKSpace.sSign(signed_cell_border) == KSpace::NEG );
            const Scalar orientation = ( flipped_border == iter_property->second.flipped ? 1 : -1 );
 
            triplets.push_back( Triplet(index_output, index_input, orientation) );
 
        }
    }
 
    typedef LinearOperator<Self, order, duality, order+1, duality> Derivative;
    Derivative _derivative(*this);
    ASSERT( _derivative.myContainer.rows() == kFormLength(order+1, duality) );
    ASSERT( _derivative.myContainer.cols() == kFormLength(order, duality) );
    _derivative.myContainer.setFromTriplets(triplets.begin(), triplets.end());
 
    if ( duality == DUAL && order*(dimEmbedded-order)%2 != 0 ) return -1 * _derivative;
    return _derivative;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Order order, DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, dimEmbedded-order, DGtal::OppositeDuality<duality>::duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::hodge() const
{
    BOOST_STATIC_ASSERT(( order >= 0 ));
    BOOST_STATIC_ASSERT(( order <= dimEmbedded ));
 
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
 
    typedef typename TLinearAlgebraBackend::Triplet Triplet;
    typedef std::vector<Triplet> Triplets;
    Triplets triplets;
 
    // iterate over output form values
    for (Index index=0; index<kFormLength(order, duality); index++)
    {
        const Cell cell = myKSpace.unsigns(myIndexSignedCells[actualOrder(order, duality)][index]);
 
        const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
        ASSERT( iter_property != myCellProperties.end() );
        ASSERT( iter_property->second.index == index );
 
        const Scalar size_ratio = ( duality == DGtal::PRIMAL ?
            iter_property->second.dual_size/iter_property->second.primal_size :
            iter_property->second.primal_size/iter_property->second.dual_size );
        triplets.push_back( Triplet(index, index, hodgeSign(cell, duality) * size_ratio) );
    }
 
    typedef LinearOperator<Self, order, duality, dimEmbedded-order, OppositeDuality<duality>::duality> Hodge;
    Hodge _hodge(*this);
    ASSERT( _hodge.myContainer.rows() == _hodge.myContainer.cols() );
    ASSERT( _hodge.myContainer.rows() == kFormLength(order, duality) );
    _hodge.myContainer.setFromTriplets(triplets.begin(), triplets.end());
 
    return _hodge;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
DGtal::VectorField<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::sharp(const DGtal::KForm<Self, 1, duality>& one_form) const
{
    ASSERT( one_form.myCalculus == this );
 
    const_cast<Self*>(this)->updateCachedOperators();
    ASSERT( !myCachedOperatorsNeedUpdate );
 
    const boost::array<SparseMatrix, dimAmbient>& sharp_operator_matrix = mySharpOperatorMatrixes[static_cast<int>(duality)];
 
    typedef VectorField<Self, duality> Field;
    Field field(*this);
    ASSERT( field.myCoordinates.cols() == dimAmbient);
    ASSERT( field.myCoordinates.rows() == kFormLength(0, duality));
    for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
        field.myCoordinates.col(direction) = sharp_operator_matrix[direction]*one_form.myContainer;
 
    return field;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality, 0, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::sharpDirectional(const DGtal::Dimension& direction) const
{
    const_cast<Self*>(this)->updateCachedOperators();
    ASSERT( !myCachedOperatorsNeedUpdate );
 
    ASSERT( direction < dimAmbient );
 
    typedef LinearOperator<Self, 1, duality, 0, duality> Operator;
    return Operator(*this, mySharpOperatorMatrixes[static_cast<int>(duality)][direction]);
}
 
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
void
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateSharpOperator()
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    ASSERT( myCachedOperatorsNeedUpdate );
 
    typedef typename TLinearAlgebraBackend::Triplet Triplet;
    typedef std::vector<Triplet> Triplets;
    typedef typename Properties::const_iterator PropertiesConstIterator;
 
    boost::array<Triplets, dimAmbient> triplets;
 
    // iterate over points
    for (Index point_index=0; point_index<kFormLength(0, duality); point_index++)
    {
        const SCell signed_point = myIndexSignedCells[actualOrder(0, duality)][point_index];
        ASSERT( myKSpace.sDim(signed_point) == actualOrder(0, duality) );
        const Scalar point_orientation = ( myKSpace.sSign(signed_point) == KSpace::POS ? 1 : -1 );
        const Cell point = myKSpace.unsigns(signed_point);
 
        typedef typename KSpace::Cells Edges;
        typedef typename Edges::const_iterator EdgesConstIterator;
        const Edges edges = ( duality == PRIMAL ? myKSpace.uUpperIncident(point) : myKSpace.uLowerIncident(point) );
        ASSERT( edges.size() <= 2*dimAmbient );
 
        // collect 1-form values over neighboring edges
        typedef boost::array< std::pair<Scalar, std::list< std::pair<Index, Scalar> > >, dimAmbient > EdgeIndexes;
        EdgeIndexes edge_indexes;
        for (EdgesConstIterator ei=edges.begin(), eie=edges.end(); ei!=eie; ei++)
        {
            const Cell edge = *ei;
            ASSERT( myKSpace.uDim(edge) == actualOrder(1, duality) );
 
            const PropertiesConstIterator edge_property_iter = myCellProperties.find(edge);
            if (edge_property_iter == myCellProperties.end())
                continue;
 
            const Index edge_index = edge_property_iter->second.index;
            const Scalar edge_length = ( duality == PRIMAL ? edge_property_iter->second.primal_size : edge_property_iter->second.dual_size );
 
            const Scalar edge_orientation = ( edge_property_iter->second.flipped ? 1 : -1 );
            const DGtal::Dimension edge_direction = edgeDirection(edge, duality); //FIXME iterate over direction
 
            edge_indexes[edge_direction].second.push_back(std::make_pair(edge_index, edge_orientation));
            edge_indexes[edge_direction].first += edge_length;
        }
 
        for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
        {
            const Scalar edge_sign = ( duality == DUAL && (direction*(dimAmbient-direction))%2 == 0 ? -1 : 1 );
            const Scalar edge_length_sum = edge_indexes[direction].first;
 
            for (typename std::list< std::pair<Index, Scalar> >::const_iterator ei=edge_indexes[direction].second.begin(), ee=edge_indexes[direction].second.end(); ei!=ee; ei++)
            {
                const Index edge_index = ei->first;
                const Scalar edge_orientation = ei->second;
                ASSERT( edge_index < static_cast<Index>(myIndexSignedCells[actualOrder(1, duality)].size()) );
                ASSERT( edge_index < kFormLength(1, duality) );
                ASSERT( edge_length_sum > 0 );
 
                triplets[direction].push_back( Triplet(point_index, edge_index, point_orientation*edge_sign*edge_orientation/edge_length_sum) );
            }
        }
    }
 
    boost::array<SparseMatrix, dimAmbient> sharp_operator_matrix;
 
    for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
    {
        sharp_operator_matrix[direction] = SparseMatrix(kFormLength(0, duality), kFormLength(1, duality));
        sharp_operator_matrix[direction].setFromTriplets(triplets[direction].begin(), triplets[direction].end());
    }
 
    mySharpOperatorMatrixes[static_cast<int>(duality)] = sharp_operator_matrix;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
DGtal::KForm<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::flat(const DGtal::VectorField<Self, duality>& vector_field) const
{
    ASSERT( vector_field.myCalculus == this );
 
    const_cast<Self*>(this)->updateCachedOperators();
    ASSERT( !myCachedOperatorsNeedUpdate );
 
    const boost::array<SparseMatrix, dimAmbient>& flat_operator_matrix = myFlatOperatorMatrixes[static_cast<int>(duality)];
 
    typedef KForm<Self, 1, duality> OneForm;
    OneForm one_form(*this);
    for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
        one_form.myContainer += flat_operator_matrix[direction]*vector_field.myCoordinates.col(direction);
 
    return one_form;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 1, duality>
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::flatDirectional(const DGtal::Dimension& direction) const
{
    const_cast<Self*>(this)->updateCachedOperators();
    ASSERT( !myCachedOperatorsNeedUpdate );
 
    ASSERT( direction < dimAmbient );
 
    typedef LinearOperator<Self, 0, duality, 1, duality> Operator;
    return Operator(*this, myFlatOperatorMatrixes[static_cast<int>(duality)][direction]);
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Duality duality>
void
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateFlatOperator()
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    ASSERT( myCachedOperatorsNeedUpdate );
 
    typedef typename TLinearAlgebraBackend::Triplet Triplet;
    typedef std::vector<Triplet> Triplets;
    typedef typename Properties::const_iterator PropertiesConstIterator;
 
    boost::array<Triplets, dimAmbient> triplets;
 
    // iterate over edges
    for (Index edge_index=0; edge_index<kFormLength(1, duality); edge_index++)
    {
        const SCell signed_edge = myIndexSignedCells[actualOrder(1, duality)][edge_index];
        ASSERT( myKSpace.sDim(signed_edge) == actualOrder(1, duality) );
        const Cell edge = myKSpace.unsigns(signed_edge);
 
        const Scalar edge_orientation = ( myKSpace.sSign(signed_edge) == KSpace::NEG ? 1 : -1 );
        const DGtal::Dimension& edge_direction = edgeDirection(edge, duality); //FIXME iterate over edge direction
        const Scalar edge_sign = ( duality == DUAL && (edge_direction*(dimAmbient-edge_direction))%2 == 0 ? -1 : 1 );
        const PropertiesConstIterator edge_property_iter = myCellProperties.find(edge);
        ASSERT( edge_property_iter != myCellProperties.end() );
        const Scalar edge_length = ( duality == PRIMAL ? edge_property_iter->second.primal_size : edge_property_iter->second.dual_size );
 
        typedef typename KSpace::Cells Points;
        const Points points = ( duality == PRIMAL ? myKSpace.uLowerIncident(edge) : myKSpace.uUpperIncident(edge) );
 
        // project vector field along edge from neighboring points
        typedef std::pair<Index, Scalar> BorderInfo;
        typedef std::list<BorderInfo> BorderInfos;
        BorderInfos border_infos;
        for (typename Points::const_iterator pi=points.begin(), pie=points.end(); pi!=pie; pi++)
        {
            const Cell point = *pi;
            ASSERT( myKSpace.uDim(point) == actualOrder(0, duality) );
 
            const PropertiesConstIterator point_property_iter = myCellProperties.find(point);
            if (point_property_iter == myCellProperties.end())
                continue;
 
            const Index point_index = point_property_iter->second.index;
            const Scalar point_orientation = ( point_property_iter->second.flipped ? -1 : 1 );
 
            border_infos.push_back(std::make_pair(point_index, point_orientation));
        }
 
        ASSERT( border_infos.size() <= 2 );
 
        for (typename BorderInfos::const_iterator bi=border_infos.begin(), bie=border_infos.end(); bi!=bie; bi++)
        {
            const Index point_index = bi->first;
            const Scalar point_orientation = bi->second;
            ASSERT( point_index < static_cast<Index>(myIndexSignedCells[actualOrder(0, duality)].size()) );
            ASSERT( point_index < kFormLength(0, duality) );
 
            triplets[edge_direction].push_back( Triplet(edge_index, point_index, point_orientation*edge_length*edge_sign*edge_orientation/border_infos.size()) );
        }
 
    }
 
    boost::array<SparseMatrix, dimAmbient> flat_operator_matrix;
 
    for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
    {
        flat_operator_matrix[direction] = SparseMatrix(kFormLength(1, duality), kFormLength(0, duality));
        flat_operator_matrix[direction].setFromTriplets(triplets[direction].begin(), triplets[direction].end());
    }
 
    myFlatOperatorMatrixes[static_cast<int>(duality)] = flat_operator_matrix;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
void
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateIndexes()
{
    if (!myIndexesNeedUpdate) return;
 
    // clear index signed cells
    for (DGtal::Dimension dim=0; dim<dimEmbedded+1; dim++)
        myIndexSignedCells[dim].clear();
 
    // compute cell index
    for (typename Properties::iterator csi=myCellProperties.begin(), csie=myCellProperties.end(); csie!=csi; csi++)
    {
        const Cell& cell = csi->first;
        const DGtal::Dimension cell_dim = myKSpace.uDim(cell);
 
        csi->second.index = myIndexSignedCells[cell_dim].size();
 
        const SCell& signed_cell = myKSpace.signs(cell, csi->second.flipped ? KSpace::NEG : KSpace::POS);
        myIndexSignedCells[cell_dim].push_back(signed_cell);
    }
 
    myIndexesNeedUpdate = false;
    myCachedOperatorsNeedUpdate = true;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
void
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateCachedOperators()
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    if (!myCachedOperatorsNeedUpdate) return;
    updateFlatOperator<PRIMAL>();
    updateFlatOperator<DUAL>();
    updateSharpOperator<PRIMAL>();
    updateSharpOperator<DUAL>();
    myCachedOperatorsNeedUpdate = false;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Properties&
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getProperties() const
{
    return myCellProperties;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
template <DGtal::Order order, DGtal::Duality duality>
const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCells&
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getIndexedSCells() const
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    return myIndexSignedCells[actualOrder(order, duality)];
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCell
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getSCell(const Order& order, const Duality& duality, const Index& index) const
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    const Order& actual_order = actualOrder(order, duality);
    const SCell& signed_cell = myIndexSignedCells[actual_order][index];
    ASSERT( myKSpace.sDim(signed_cell) == actual_order );
    return signed_cell;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
bool
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::containsCell(const Cell& cell) const
{
    return myCellProperties.find(cell) != myCellProperties.end();
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
bool
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::isCellFlipped(const Cell& cell) const
{
    const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
    ASSERT( iter_property != myCellProperties.end() );
    return iter_property->second.flipped;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getCellIndex(const Cell& cell) const
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
    ASSERT( iter_property != myCellProperties.end() );
    return iter_property->second.index;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin() const
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    return myCellProperties.begin();
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end() const
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    return myCellProperties.end();
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin()
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    return myCellProperties.begin();
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end()
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    return myCellProperties.end();
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::kFormLength(const DGtal::Order& order, const DGtal::Duality& duality) const
{
    ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
    return myIndexSignedCells[actualOrder(order, duality)].size();
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
DGtal::Order
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::actualOrder(const DGtal::Order& order, const DGtal::Duality& duality) const
{
    return duality == PRIMAL ? order : dimEmbedded-order;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Scalar
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::hodgeSign(const Cell& cell, const DGtal::Duality& duality) const
{
    if (duality == PRIMAL) return 1;
 
    const DGtal::Dimension& primal_dim = myKSpace.uDim(cell);
    return (dimEmbedded-primal_dim)*primal_dim % 2 != 0 ? -1 : 1;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
typename DGtal::Dimension
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::edgeDirection(const Cell& cell, const DGtal::Duality& duality) const
{
    ASSERT( myKSpace.uDim(cell) == actualOrder(1, duality) );
 
    typename KSpace::DirIterator di = myKSpace.uDirs(cell);
    if (duality == DUAL) di = myKSpace.uOrthDirs(cell);
 
    ASSERT( di != 0 );
    const DGtal::Dimension direction = *di;
    ++di;
    //ASSERT( !(di != 0) ); //FIXME multiple edge direction with embedding
    ASSERT( direction < dimAmbient );
 
    return direction;
}
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
bool
DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::isValid() const
{
    return true;
}
 
///////////////////////////////////////////////////////////////////////////////
// Implementation of inline functions                                        //
 
template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
std::ostream&
DGtal::operator<<(std::ostream & out, const DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>& object)
{
  object.selfDisplay(out);
  return out;
}
 
//                                                                           //
///////////////////////////////////////////////////////////////////////////////