derivatives.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */



#ifndef SLIP_DERIVATIVES_HPP
#define SLIP_DERIVATIVES_HPP
#include <iostream>
#include <cassert>
#include <vector>
#include "linear_algebra.hpp"
#include "linear_algebra_traits.hpp"
#include "convolution.hpp"
#include "Vector.hpp"
#include "Matrix.hpp"
#include "Point2d.hpp"
#include "Point3d.hpp"
#include "Point4d.hpp"
#include "neighbors.hpp"

namespace slip
{

enum SPATIAL_DIRECTION
{
   X_DIRECTION, 
   Y_DIRECTION, 
   Z_DIRECTION
};

  /* @{ */

  template<class Container2D>
  inline
  void finite_diff_coef(const std::size_t sch_ord, 
                                                const std::size_t sch_shift, 
                                                Container2D& Dinv)
  {
    //Precondition
  assert(sch_ord >= sch_shift);
  // Internal Variables
  slip::Vector<typename Container2D::value_type> tmp(sch_ord+1);
  Container2D D(sch_ord+1,sch_ord+1,0.0);
     
  for(std::size_t i = 0; i <= sch_ord; ++i)
  {
    D[i][0] = 1.0;
    tmp[i] = typename Container2D::value_type(int(i) - int(sch_shift));
  }
  for(std::size_t i = 0; i <= sch_ord; ++i)
    {
      for(std::size_t j = 1; j <= sch_ord; ++j)
                {
                  D[i][j] = D[i][j-1] * tmp[i] / typename Container2D::value_type(j);
                }
    }
  // Result container
  Dinv.resize(sch_ord+1,sch_ord+1); //finite differences coefficient matrix
  //matrix inversion
  slip::inverse(D,Dinv);
}

  template<class Container2D>
  inline
  void finite_diff_coef(const std::size_t sch_ord, 
                                                const std::size_t sch_shift, 
                                                const typename Container2D::value_type& h,
                                                Container2D& Dinv)
  {
    //Precondition
  assert(sch_ord >= sch_shift);
  assert(h != typename Container2D::value_type(0));
  typedef typename Container2D::value_type value_type;
  // Internal Variables
  slip::Vector<value_type> tmp(sch_ord+1);
  Container2D D(sch_ord+1,sch_ord+1,0.0);
     
  for(std::size_t i = 0; i <= sch_ord; ++i)
  {
    D[i][0] = 1.0;
    tmp[i] = value_type(int(i) - int(sch_shift));
  }
  for(std::size_t i = 0; i <= sch_ord; ++i)
    {
      for(std::size_t j = 1; j <= sch_ord; ++j)
                {
                  D[i][j] = D[i][j-1] * tmp[i] / value_type(j);
                }
    }
  // Result container
  Dinv.resize(sch_ord+1,sch_ord+1); //finite differences coefficient matrix
  //matrix inversion
  slip::inverse(D,Dinv);
  value_type inv_h = value_type(1)/h;
  for(std::size_t i = 1; i <= sch_ord; ++i)
    {
      for(std::size_t j = 0; j <= sch_ord; ++j)
                {
                  Dinv[i][j] *= inv_h;
                }
      inv_h*=inv_h;
    }
  
}

  template<typename KernelsIterator>
  inline
  void finite_diff_kernels(const std::size_t der_ord,
                                                   const std::size_t sch_ord, 
                                                   const std::size_t sch_shift, 
                                                   std::vector<KernelsIterator>& kernels_first)
  {
    //central kernel
    slip::Matrix<typename std::iterator_traits<KernelsIterator>::value_type> Dinv(sch_ord+1,sch_ord+1,0.0);
    slip::finite_diff_coef(sch_ord,sch_shift,Dinv);
    std::copy(Dinv.row_rbegin(der_ord),Dinv.row_rend(der_ord),kernels_first[0]);
    //first border kernels
     for (std::size_t bord_shift = 0; bord_shift < sch_shift; ++bord_shift)
                {
                  //generate finite differences coefficients
                  slip::finite_diff_coef(sch_ord,bord_shift,Dinv);
                  //extract derivation kernel
                  std::copy(Dinv.row_rbegin(der_ord),Dinv.row_rend(der_ord),kernels_first[bord_shift + 1]);

                }
     //last border kernels
     for (std::size_t bord_shift = sch_shift + 1; bord_shift <= sch_ord; ++bord_shift)
       {
                 //generate finite differences coefficients
                 slip::finite_diff_coef(sch_ord,bord_shift,Dinv);
                 //extract derivation kernel
                 std::copy(Dinv.row_rbegin(der_ord),Dinv.row_rend(der_ord),kernels_first[bord_shift]);
       }
  }
  
  template<typename KernelsIterator>
  inline
  void finite_diff_kernels(const std::size_t der_ord,
                                                   const std::size_t sch_ord, 
                                                   const std::size_t sch_shift,
                                                   const typename std::iterator_traits<KernelsIterator>::value_type& h,
                                                   std::vector<KernelsIterator>& kernels_first)
  {
    //central kernel
    slip::Matrix<typename std::iterator_traits<KernelsIterator>::value_type> Dinv(sch_ord+1,sch_ord+1,0.0);
    slip::finite_diff_coef(sch_ord,sch_shift,h,Dinv);
    std::copy(Dinv.row_rbegin(der_ord),Dinv.row_rend(der_ord),kernels_first[0]);
    //first border kernels
     for (std::size_t bord_shift = 0; bord_shift < sch_shift; ++bord_shift)
                {
                  //generate finite differences coefficients
                  slip::finite_diff_coef(sch_ord,bord_shift,h,Dinv);
                  //extract derivation kernel
                  std::copy(Dinv.row_rbegin(der_ord),Dinv.row_rend(der_ord),kernels_first[bord_shift + 1]);

                }
     //last border kernels
     for (std::size_t bord_shift = sch_shift + 1; bord_shift <= sch_ord; ++bord_shift)
       {
                 //generate finite differences coefficients
                 slip::finite_diff_coef(sch_ord,bord_shift,h,Dinv);
                 //extract derivation kernel
                 std::copy(Dinv.row_rbegin(der_ord),Dinv.row_rend(der_ord),kernels_first[bord_shift]);
       }
  }

  template<typename SrcIter, typename Container2D, typename ResIter>
  inline
  void derivative(SrcIter first, SrcIter last, 
                                  const std::size_t der_ord, 
                                  const std::size_t sch_ord, 
                                  const std::size_t sch_shift,
                                  const Container2D& M,
                                  ResIter result)
  {
    assert(sch_ord >= der_ord);
    assert(sch_shift <= sch_ord);
    //- derivation kernel
    slip::Vector<double> K(sch_ord+1);
    //extract derivation kernel (reversed for convolution algorithm)
    for(std::size_t i=0; i<=sch_ord; ++i)
                K[sch_ord-i]=M[der_ord][i];

    slip::same_convolution(first,last,
                                                   K.begin(),K.end(),sch_ord-sch_shift,sch_shift,
                                                   result,
                                                   slip::BORDER_TREATMENT_AVOID);
    // Border Treatment                                 //
    // A variable derivation kernel is applied :        //
    // -scheme shift varies from 0 to sch_shift         //
    //  on the first elements                           //
    // -scheme shift varies from sch_shift+1 to sch_ord //
    //  on the last elements                            //

    // Internal Variables
     //- derivation kernel
    slip::Matrix<double> Dinv(sch_ord+1,sch_ord+1);
    // beginning of block : shift from 0 to sch_shift-1
      for (std::size_t bord_shift = 0; bord_shift < sch_shift; ++bord_shift)
                {
        //generate finite differences coefficients
                  slip::finite_diff_coef(sch_ord,bord_shift,Dinv);
        //extract derivation kernel (reversed for convolution algorithm)
                  for(std::size_t i = 0; i <= sch_ord; ++i)
                    K[sch_ord - i] = Dinv[der_ord][i];
                  //compute the convolution on the left border
                  slip::valid_convolution(first,first+(sch_ord + 1),
                                                                  K.begin(),K.end(),(sch_ord - bord_shift), 
                                                                  bord_shift,
                                                                  result + bord_shift);
                }  

    // end of block : shift from sch_shift+1 to sch_ord
      for (std::size_t bord_shift = sch_shift + 1; bord_shift <= sch_ord; ++bord_shift)
                {
                  //generate finite differences coefficients
                  slip::finite_diff_coef(sch_ord,bord_shift,Dinv);
                  //extract derivation kernel (reversed for convolution algorithm)
                  for(std::size_t i = 0; i <= sch_ord; ++i)
                    K[sch_ord - i] = Dinv[der_ord][i];

                  //compute the convolution at the end of block
                  slip::valid_convolution(last - (1 + sch_ord),
                                                                  last,
                                                                  K.begin(),K.end(),(sch_ord - bord_shift),bord_shift,
                                                                  (result + (last-first))-(1 + sch_ord) + bord_shift); 
                }
  }

  template<typename SrcIter, typename KernelsIterator, typename ResIter>
  inline
  void derivative(SrcIter first, SrcIter last, 
                                  const std::size_t der_ord, 
                                  const std::size_t sch_ord, 
                                  const std::size_t sch_shift,
                                  const std::vector<KernelsIterator>& kernels_first,
                                  ResIter result)
  {
    assert(sch_ord >= der_ord);
    assert(sch_shift <= sch_ord);

    std::size_t sch_ord_plus_one = sch_ord + 1;
    slip::same_convolution(first,last,
                                                   kernels_first[0],kernels_first[0] + sch_ord_plus_one,
                                                   sch_ord-sch_shift,
                                                   sch_shift,
                                                   result,
                                                   slip::BORDER_TREATMENT_AVOID);
    // Border Treatment                                 //
    // A variable derivation kernel is applied :        //
    // -scheme shift varies from 0 to sch_shift         //
    //  on the first elements                           //
    // -scheme shift varies from sch_shift+1 to sch_ord //
    //  on the last elements                            //

    // beginning of block : shift from 0 to sch_shift-1
    for (std::size_t bord_shift = 0; bord_shift < sch_shift; ++bord_shift)
      {
                //compute the convolution on the left border
                slip::valid_convolution(first,first + sch_ord_plus_one,
                                                                kernels_first[bord_shift+1],
                                                                kernels_first[bord_shift+1]+ sch_ord_plus_one,
                                                                (sch_ord - bord_shift), 
                                                                bord_shift,
                                                                result + bord_shift);
      }  
    
    // end of block : shift from sch_shift+1 to sch_ord
    for (std::size_t bord_shift = sch_shift + 1; bord_shift <= sch_ord; ++bord_shift)
      {
                //compute the convolution at the end of block
                slip::valid_convolution(last - sch_ord_plus_one,
                                                                last,
                                                                kernels_first[bord_shift],
                                                                kernels_first[bord_shift]+ sch_ord_plus_one,
                                                                (sch_ord - bord_shift),
                                                                bord_shift,
                                                                (result + (last-first))-sch_ord_plus_one + bord_shift); 
                }
  }

  
  template<typename Iterator2d1, typename Iterator2d2, typename T>
  inline
  void derivative_2d(Iterator2d1 I_up,
                                     Iterator2d1 I_bot,
                                     const slip::Point2d<T>& grid_step,
                                     slip::SPATIAL_DIRECTION der_dir,
                                     const std::size_t der_order, 
                                     const std::size_t sch_order,
                                     Iterator2d2 R_up,
                                     Iterator2d2 R_bot) 
  {
    assert(sch_order >= der_order);
    assert(   ((der_dir == X_DIRECTION) && (int((I_bot-I_up)[0]) > int(sch_order)))
                      ||((der_dir == Y_DIRECTION) && (int((I_bot-I_up)[1]) > int(sch_order))) );
    assert((sch_order / 2) <= sch_order);
    assert((R_bot - R_up) == (I_bot-I_up));
    typedef typename Iterator2d1::value_type value_type;
    typedef typename Iterator2d1::size_type size_type;
    size_type sch_shift = sch_order / 2;
   
    
    //computes all kernels
    slip::Matrix<T> kernels(sch_order + 1, sch_order + 1);
    std::vector<typename slip::Matrix<T>::iterator> kernels_iterators(sch_order + 1);   
    for(size_type i = 0; i < (sch_order + 1); ++i)
      {
                kernels_iterators[i] = kernels.row_begin(i);
      }
    
    if(der_dir == X_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                      grid_step[0],
                                                      kernels_iterators);
                size_type rows = size_type((I_bot-I_up)[0]);
                for(size_type i = 0; i < rows; ++i)
                  {
                    slip::derivative(I_up.row_begin(i),
                                                     I_up.row_end(i), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernels_iterators,
                                                     R_up.row_begin(i));
                  }
      }
    
    if(der_dir == Y_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                      grid_step[1],
                                                      kernels_iterators);
                size_type cols = size_type((I_bot-I_up)[1]);
                for(size_type j = 0; j < cols; ++j)
                  {
                    slip::derivative(I_up.col_rbegin(j),
                                                     I_up.col_rend(j), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernels_iterators,
                                                     R_up.col_rbegin(j));
                    
                  }
      }
  }

  template<typename RandomAccessIterator3d1, 
                   typename RandomAccessIterator3d2, 
                   typename T>
  inline
  void derivative_3d(RandomAccessIterator3d1 I_fup,
                                     RandomAccessIterator3d1 I_back_bot,
                                     const slip::Point3d<T>& grid_step,
                                     slip::SPATIAL_DIRECTION der_dir,
                                     const std::size_t der_order, 
                                     const std::size_t sch_order,
                                     RandomAccessIterator3d2 R_fup,
                                     RandomAccessIterator3d2 R_back_bot) 
  {
    assert(sch_order >= der_order);
    assert(   ((der_dir == X_DIRECTION) && (int((I_back_bot-I_fup)[2]) > int(sch_order)))
                      ||((der_dir == Y_DIRECTION) && (int((I_back_bot-I_fup)[1]) > int(sch_order))) 
                      ||((der_dir == Z_DIRECTION) && (int((I_back_bot-I_fup)[0]) > int(sch_order))));
    assert((sch_order / 2) <= sch_order);
    assert((R_back_bot - R_fup) == (I_back_bot-I_fup));
    typedef typename RandomAccessIterator3d1::value_type value_type;
    typedef typename RandomAccessIterator3d1::size_type size_type;
    size_type sch_shift = sch_order / 2;
   
    
    //computes all kernels
    slip::Matrix<T> kernels(sch_order + 1, sch_order + 1);
    std::vector<typename slip::Matrix<T>::iterator> kernels_iterators(sch_order + 1);
    std::vector<typename slip::Matrix<T>::reverse_iterator> rkernels_iterators(sch_order + 1);
    for(size_type i = 0; i < (sch_order + 1); ++i)
      {
                kernels_iterators[i] = kernels.row_begin(i);
                rkernels_iterators[i] = kernels.row_rbegin(i);
      }
    
    if(der_dir == X_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step[0],
                                                                  kernels_iterators);
                const size_type rows   = static_cast<std::size_t>((I_back_bot-I_fup)[1]);
                const size_type slices = static_cast<std::size_t>((I_back_bot-I_fup)[0]);
                for(size_type k = 0; k < slices; ++k)
                  {
                    for(size_type i = 0; i < rows; ++i)
                      {
                                slip::derivative(I_fup.row_begin(k,i),
                                                                 I_fup.row_end(k,i), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 kernels_iterators,
                                                                 R_fup.row_begin(k,i));
                      }
                  }
      }
    
    if(der_dir == Y_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step[1],
                                                                  kernels_iterators);
                const size_type cols   = static_cast<std::size_t>((I_back_bot-I_fup)[2]);
                const size_type slices = static_cast<std::size_t>((I_back_bot-I_fup)[0]);

                for(size_type k = 0; k < slices; ++k)
                  {
                    for(size_type j = 0; j < cols; ++j)
                      {
                                slip::derivative(I_fup.col_begin(k,j),
                                                                 I_fup.col_end(k,j), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 rkernels_iterators,
                                                                 R_fup.col_begin(k,j));
                      }
                  }
      }
    if(der_dir == Z_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step[2],
                                                                  kernels_iterators);
                const size_type rows   = static_cast<std::size_t>((I_back_bot-I_fup)[1]);
                const size_type cols   = static_cast<std::size_t>((I_back_bot-I_fup)[2]);

                for(size_type i = 0; i < rows; ++i)
                  {
                    for(size_type j = 0; j < cols; ++j)
                      {
                                slip::derivative(I_fup.slice_begin(i,j),
                                                                 I_fup.slice_end(i,j), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 rkernels_iterators,
                                                                 R_fup.slice_begin(i,j));
                      }
                  }
      }
  }


enum FOUR_DIM_DIRECTION
{
                T_DIRECTION,
                S_DIRECTION,
                R_DIRECTION,
                C_DIRECTION
};

/* @{ */

template<typename RandomAccessIterator4d1,
typename RandomAccessIterator4d2,
typename T>
inline
void derivative_4d(RandomAccessIterator4d1 I_ffup,
                                RandomAccessIterator4d1 I_last_back_bot,
                                const slip::Point4d<T>& grid_step,
                                slip::FOUR_DIM_DIRECTION der_dir,
                                const std::size_t der_order,
                                const std::size_t sch_order,
                                RandomAccessIterator4d2 R_ffup,
                                RandomAccessIterator4d2 R_last_back_bot)
{
                assert(sch_order >= der_order);
                assert(   ((der_dir == T_DIRECTION) && (int((I_last_back_bot-I_ffup)[0]) > int(sch_order)))
                                                ||((der_dir == S_DIRECTION) && (int((I_last_back_bot-I_ffup)[3]) > int(sch_order)))
                                                ||((der_dir == R_DIRECTION) && (int((I_last_back_bot-I_ffup)[2]) > int(sch_order)))
                                                ||((der_dir == C_DIRECTION) && (int((I_last_back_bot-I_ffup)[1]) > int(sch_order))));
                assert((sch_order / 2) <= sch_order);
                assert((R_last_back_bot - R_ffup) == (I_last_back_bot-I_ffup));
                typedef typename RandomAccessIterator4d1::value_type value_type;
                typedef typename RandomAccessIterator4d1::size_type size_type;
                size_type sch_shift = sch_order / 2;


                //computes all kernels
                slip::Matrix<T> kernels(sch_order + 1, sch_order + 1);
                std::vector<typename slip::Matrix<T>::iterator> kernels_iterators(sch_order + 1);
                std::vector<typename slip::Matrix<T>::reverse_iterator> rkernels_iterators(sch_order + 1);
                for(size_type i = 0; i < (sch_order + 1); ++i)
                {
                                kernels_iterators[i] = kernels.row_begin(i);
                                rkernels_iterators[i] = kernels.row_rbegin(i);
                }

                if(der_dir == T_DIRECTION)
                {
                                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                grid_step[0],
                                                                kernels_iterators);
                                const size_type rows   = static_cast<std::size_t>((I_last_back_bot-I_ffup)[2]);
                                const size_type cols   = static_cast<std::size_t>((I_last_back_bot-I_ffup)[3]);
                                const size_type slices = static_cast<std::size_t>((I_last_back_bot-I_ffup)[1]);
                                for(size_type k = 0; k < slices; ++k)
                                {
                                                for(size_type i = 0; i < rows; ++i)
                                                {
                                                                for(size_type j = 0; j < cols; ++j)
                                                                {
                                                                                slip::derivative(I_ffup.slab_begin(k,i,j),
                                                                                                                I_ffup.slab_end(k,i,j),
                                                                                                                der_order,
                                                                                                                sch_order,
                                                                                                                sch_shift,
                                                                                                                kernels_iterators,
                                                                                                                R_ffup.slab_begin(k,i,j));
                                                                }
                                                }
                                }
                }

                if(der_dir == S_DIRECTION)
                {
                                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                grid_step[1],
                                                                kernels_iterators);
                                const size_type rows   = static_cast<std::size_t>((I_last_back_bot-I_ffup)[2]);
                                const size_type slices = static_cast<std::size_t>((I_last_back_bot-I_ffup)[1]);
                                const size_type slabs = static_cast<std::size_t>((I_last_back_bot-I_ffup)[0]);
                                for(size_type t = 0; t < slabs; ++t)
                                {
                                                for(size_type k = 0; k < slices; ++k)
                                                {
                                                                for(size_type i = 0; i < rows; ++i)
                                                                {
                                                                                slip::derivative(I_ffup.row_begin(t,k,i),
                                                                                                                I_ffup.row_end(t,k,i),
                                                                                                                der_order,
                                                                                                                sch_order,
                                                                                                                sch_shift,
                                                                                                                kernels_iterators,
                                                                                                                R_ffup.row_begin(t,k,i));
                                                                }
                                                }
                                }
                }

                if(der_dir == R_DIRECTION)
                {
                                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                grid_step[2],
                                                                kernels_iterators);
                                const size_type cols   = static_cast<std::size_t>((I_last_back_bot-I_ffup)[3]);
                                const size_type slices = static_cast<std::size_t>((I_last_back_bot-I_ffup)[1]);
                                const size_type slabs = static_cast<std::size_t>((I_last_back_bot-I_ffup)[0]);
                                for(size_type t = 0; t < slabs; ++t)
                                {
                                                for(size_type k = 0; k < slices; ++k)
                                                {
                                                                for(size_type j = 0; j < cols; ++j)
                                                                {
                                                                                slip::derivative(I_ffup.col_begin(t,k,j),
                                                                                                                I_ffup.col_end(t,k,j),
                                                                                                                der_order,
                                                                                                                sch_order,
                                                                                                                sch_shift,
                                                                                                                rkernels_iterators,
                                                                                                                R_ffup.col_begin(t,k,j));
                                                                }
                                                }
                                }
                }
                if(der_dir == C_DIRECTION)
                {
                                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                grid_step[3],
                                                                kernels_iterators);
                                const size_type rows   = static_cast<std::size_t>((I_last_back_bot-I_ffup)[2]);
                                const size_type cols   = static_cast<std::size_t>((I_last_back_bot-I_ffup)[3]);
                                const size_type slabs = static_cast<std::size_t>((I_last_back_bot-I_ffup)[0]);
                                for(size_type t = 0; t < slabs; ++t)
                                {
                                                for(size_type i = 0; i < rows; ++i)
                                                {
                                                                for(size_type j = 0; j < cols; ++j)
                                                                {
                                                                                slip::derivative(I_ffup.slice_begin(t,i,j),
                                                                                                                I_ffup.slice_end(t,i,j),
                                                                                                                der_order,
                                                                                                                sch_order,
                                                                                                                sch_shift,
                                                                                                                rkernels_iterators,
                                                                                                                R_ffup.slice_begin(t,i,j));
                                                                }
                                                }
                                }
                }
}




/* @} */


  //----------------------------------------------------------
  // local derivatives
  //----------------------------------------------------------
  
template<typename Real> 
inline
Real minmod(const Real gplus, const Real gminus)
{
  Real m = Real();

  if (gplus*gminus > Real())
   {
     Real absgplus  = std::abs(gplus);
     Real absgminus = std::abs(gminus);
     m = std::min(absgplus,absgminus);
     m = m * m;
   }
 else
   {
     m = Real(); 
   }
  return m;
}

template<typename Real>
inline
Real max_brockett_maragos_plus(const Real gplus, const Real gminus)
{

Real max = std::max(gplus,-gminus);
max = std::max(max,Real(0.0));
return max * max;
}

template<typename Real>
inline
Real max_brockett_maragos_minus(const Real gplus, const Real gminus)
{

  Real max = std::max(Real(0.0),gminus);
  max = std::max(max,gplus);
  return max * max;
}


template<typename Real>
inline
Real osher_sethian_plus(const Real gplus, const Real gminus)
{
  Real  max = std::max(gplus,Real(0.0));
  Real  min = std::min(gminus,Real(0.0));

  max = max * max;
  min = min * min;

  return min + max;
}

template<typename Real>
inline
Real osher_sethian_minus(const Real gplus, const Real gminus)
{


  Real max = std::max(gminus,Real(0.0));
  Real min = std::min(gplus,Real(0.0));
  
  max  = max * max;
  min  = min * min;

  return min + max;
}

//
// dx local derivatives
//
template<typename Real,typename>
struct __dx
{
  template <typename _II>
  static Real
  dx(_II it)
  {
    Real left    = *(it - 1);
    Real right   = *(it + 1);
    
    return Real(0.5) * (right - left);
  }
};

template<typename Real>
struct __dx<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dx(_II it)
  {
    Real left    = *(it - 1);
    Real right   = *(it + 1);
    
    return Real(0.5) * (right - left);
  }
};

template<typename Real>
struct __dx<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dx(_II it)
  {
    Real left    = it[slip::n_8c[4]];
    Real right   = it[slip::n_8c[0]];

  return Real(0.5) * (right - left);
  }
};

template<typename Real>
struct __dx<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dx(_II it)
  {
    Real left    = it[slip::n_6c[3]];
    Real right   = it[slip::n_6c[1]];

  return Real(0.5) * (right - left);
  }
};

template<typename Real>
struct __dx<Real,std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dx(_II it)
                {
                                Real left    = it[slip::n_4d_8c[4]];
                                Real right   = it[slip::n_4d_8c[2]];

                                return Real(0.5) * (right - left);
                }
};


template<typename Real> struct Dx 
{
  
  Dx(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dx<Real,_Category>::dx(it);
  }
};



//
// dy local derivatives
//
template<typename Real,typename>
struct __dy
{
  template <typename _II>
  static Real
  dy(_II it)
  {
    Real left    = *(it - 1);
    Real right   = *(it + 1);
    
    return Real(0.5) * (right - left);
  }
};

template<typename Real>
struct __dy<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dy(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __dy<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dy(_II it)
  {
    Real top      = it[slip::n_8c[2]];
    Real bottom   = it[slip::n_8c[6]];

  return Real(0.5) * (top - bottom);
  }
};

template<typename Real>
struct __dy<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dy(_II it)
  {
    Real top    = it[slip::n_6c[2]];
    Real bottom = it[slip::n_6c[4]];

  return Real(0.5) * (top - bottom);
  }
};

template<typename Real>
struct __dy<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dy(_II it)
                {
                                Real top    = it[slip::n_4d_8c[3]];
                                Real bottom = it[slip::n_4d_8c[5]];

                                return Real(0.5) * (top - bottom);
                }
};

template<typename Real> struct Dy
{
 
  Dy(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dy<Real,_Category>::dy(it);
  }
};


//
// dz local derivatives
//
template<typename Real,typename>
struct __dz
{
  template <typename _II>
  static Real
  dz(_II it)
  {
    Real left    = *(it - 1);
    Real right   = *(it + 1);
    
    return Real(0.5) * (right - left);
  }
};

template<typename Real>
struct __dz<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dz(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __dz<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dz(_II it)
  {
    Real top      = it[slip::n_8c[5]];
    Real bottom   = it[slip::n_8c[0]];

  return Real(0.5) * (top - bottom);
  }
};

template<typename Real>
struct __dz<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dz(_II it)
  {
    Real top    = it[slip::n_6c[0]];
    Real bottom = it[slip::n_6c[5]];

  return Real(0.5) * (top - bottom);
  }
};

template<typename Real>
struct __dz<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dz(_II it)
                {
                                Real top    = it[slip::n_4d_8c[1]];
                                Real bottom = it[slip::n_4d_8c[6]];

                                return Real(0.5) * (top - bottom);
                }
};


template<typename Real> struct Dz
{
  
  Dz(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dz<Real,_Category>::dz(it);
  }
};


template<typename Real,typename>
struct __dt
{
                template <typename _II>
                static Real
                dt(_II it)
                {
                                Real prev = *(it - 1);
                                Real next = *(it + 1);

                                return Real(0.5) * (next - prev);
                }
};

template<typename Real>
struct __dt<Real,std::random_access_iterator_tag>
{
                template <typename _II>
                static Real
                dt(_II it)
                {
                                return __dx<Real,std::random_access_iterator_tag>::dx(it);
                }
};

template<typename Real>
struct __dt<Real, std::random_access_iterator2d_tag>
{
                template <typename _II>
                static Real
                dt(_II it)
                {
                                Real prev = it[slip::n_8c[5]];
                                Real next = it[slip::n_8c[0]];

                                return Real(next - prev);
                }
};

template<typename Real>
struct __dt<Real, std::random_access_iterator3d_tag>
{
                template <typename _II>
                static Real
                dt(_II it)
                {
                                Real prev    = it[slip::n_26c[5]];
                                Real next = it[slip::n_26c[19]];

                                return Real(next - prev);
                }
};

template<typename Real>
struct __dt<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dt(_II it)
                {
                                Real prev = it[slip::n_4d_8c[0]];
                                Real next = it[slip::n_4d_8c[7]];

                                return Real(0.5) * (next - prev);
                }
};

template<typename Real> struct Dt
{

                Dt(){}
                template<typename Iterator>
                Real operator() (Iterator& it) const
                {
                                typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
                                return __dt<Real,_Category>::dt(it);
                }
};



//
// forward dx
//
template<typename Real,typename>
struct __dx_fwd
{
  template <typename _II>
  static Real
  dx_fwd(_II it)
  {
    Real center  = *it;
    Real right   = *(it + 1);
    
    return (right - center);
  }
};

template<typename Real>
struct __dx_fwd<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dx_fwd(_II it)
  {
    Real center  = *it;
    Real right   = *(it + 1);
    
    return (right - center);
   
  }
};

template<typename Real>
struct __dx_fwd<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dx_fwd(_II it)
  {
    Real center  = *it;
    Real right   = it[slip::n_8c[0]];

  return (right - center);
  }
};

template<typename Real>
struct __dx_fwd<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dx_fwd(_II it)
  {
    Real center  = *it;
    Real right   = it[slip::n_6c[1]];

  return (right - center);
  }
};

template<typename Real>
struct __dx_fwd<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dx_fwd(_II it)
                {
                                Real center  = *it;
                                Real right   = it[slip::n_4d_8c[2]];

                                return (right - center);
                }
};


template<typename Real> struct DxFwd
{
  
  DxFwd(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dx_fwd<Real,_Category>::dx_fwd(it);
  }
};

//
// forward dy
//
template<typename Real,typename>
struct __dy_fwd
{
  template <typename _II>
  static Real
  dy_fwd(_II it)
  {
    Real center  = *it;
    Real right   = *(it + 1);
    
    return (right - center);
  }
};

template<typename Real>
struct __dy_fwd<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dy_fwd(_II it)
  {
    return __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
  }
};

template<typename Real>
struct __dy_fwd<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dy_fwd(_II it)
  {
    Real top      = it[slip::n_8c[2]];
    Real center   = *it;

  return (top - center);
  }
};

template<typename Real>
struct __dy_fwd<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dy_fwd(_II it)
  {
    Real top    = it[slip::n_6c[2]];
    Real center = *it;

  return (top - center);
  }
};

template<typename Real>
struct __dy_fwd<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dy_fwd(_II it)
                {
                                Real top    = it[slip::n_4d_8c[3]];
                                Real center = *it;

                                return (top - center);
                }
};


template<typename Real> struct DyFwd
{
  
  DyFwd(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dy_fwd<Real,_Category>::dy_fwd(it);
  }
};

//
// dz forward
//
template<typename Real,typename>
struct __dz_fwd
{
  template <typename _II>
  static Real
  dz_fwd(_II it)
  {
    Real center  = *it;
    Real right   = *(it + 1);
    
    return (right - center);
  }
};

template<typename Real>
struct __dz_fwd<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dz_fwd(_II it)
  {
    return __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
  }
};

template<typename Real>
struct __dz_fwd<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dz_fwd(_II it)
  {
    Real top      = it[slip::n_8c[5]];
    Real center   = *it;

  return (top - center);
  }
};

template<typename Real>
struct __dz_fwd<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dz_fwd(_II it)
  {
    Real top    = it[slip::n_6c[0]];
    Real center = *it;

  return (top - center);
  }
};

template<typename Real>
struct __dz_fwd<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dz_fwd(_II it)
                {
                                Real top    = it[slip::n_4d_8c[1]];
                                Real center = *it;

                                return (top - center);
                }
};

template<typename Real> struct DzFwd 
{
 
  DzFwd(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dz_fwd<Real,_Category>::dz_fwd(it);
  }
};


template<typename Real,typename>
struct __dt_fwd
{
                template <typename _II>
                static Real
                dt_fwd(_II it)
                {
                                Real center  = *it;
                                Real next   = *(it + 1);

                                return (next - center);
                }
};

template<typename Real>
struct __dt_fwd<Real,std::random_access_iterator_tag>
{
                template <typename _II>
                static Real
                dt_fwd(_II it)
                {
                                return __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
                }
};

template<typename Real>
struct __dt_fwd<Real, std::random_access_iterator2d_tag>
{
                template <typename _II>
                static Real
                dt_fwd(_II it)
                {
                                Real center = *it;
                                Real next = it[slip::n_8c[0]];

                                return Real (next - center);
                }
};

template<typename Real>
struct __dt_fwd<Real, std::random_access_iterator3d_tag>
{
                template <typename _II>
                static Real
                dt_fwd(_II it)
                {
                                Real center    = *it;
                                Real next = it[slip::n_26c[19]];

                                return Real (next - center);
                }
};

template<typename Real>
struct __dt_fwd<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dt_fwd(_II it)
                {
                                Real next    = it[slip::n_4d_8c[7]];
                                Real center = *it;

                                return (next - center);
                }
};

template<typename Real> struct DtFwd
{

                DtFwd(){}
                template<typename Iterator>
                Real operator() (Iterator& it) const
                {
                                typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
                                return __dt_fwd<Real,_Category>::dt_fwd(it);
                }
};

//
// backward dx
//
template<typename Real,typename>
struct __dx_bck
{
  template <typename _II>
  static Real
  dx_bck(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    
    return (center - left);
  }
};

template<typename Real>
struct __dx_bck<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dx_bck(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    
    return (center - left);
  }
};

template<typename Real>
struct __dx_bck<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dx_bck(_II it)
  {
    Real left    = it[slip::n_8c[4]];
    Real center  = *it;
    
    return (center - left);
  }
};

template<typename Real>
struct __dx_bck<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dx_bck(_II it)
  {
    Real left    = it[slip::n_6c[3]];
    Real center  = *it;
    
    return (center - left);
  }
};


template<typename Real>
struct __dx_bck<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dx_bck(_II it)
                {
                                Real left    = it[slip::n_4d_8c[4]];
                                Real center  = *it;

                                return (center - left);
                }
};


template<typename Real> struct DxBck
{
 
  DxBck(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dx_bck<Real,_Category>::dx_bck(it);
  }
};

//
// Backward dy
//
template<typename Real,typename>
struct __dy_bck
{
  template <typename _II>
  static Real
  dy_bck(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    
    return (center - left);
  }
};

template<typename Real>
struct __dy_bck<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dy_bck(_II it)
  {
    return __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
  }
};

template<typename Real>
struct __dy_bck<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dy_bck(_II it)
  {
    Real center   = *it;
    Real bottom   = it[slip::n_8c[6]];

  return (center - bottom);
  }
};

template<typename Real>
struct __dy_bck<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dy_bck(_II it)
  {
    Real center = *it;
    Real bottom = it[slip::n_6c[4]];

  return Real(0.5) * (center - bottom);
  }
};

template<typename Real>
struct __dy_bck<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dy_bck(_II it)
                {
                                Real center = *it;
                                Real bottom = it[slip::n_4d_8c[5]];

                                return (center - bottom);
                }
};


template<typename Real> struct DyBck
{
  
  DyBck(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dy_bck<Real,_Category>::dy_bck(it);
  }
};

//
// Backward dz
//

template<typename Real,typename>
struct __dz_bck
{
  template <typename _II>
  static Real
  dz_bck(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    
    return (left - center);
  }
};

template<typename Real>
struct __dz_bck<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dz_bck(_II it)
  {
    return __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
  }
};

template<typename Real>
struct __dz_bck<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dz_bck(_II it)
  {
    Real center   = *it;
    Real bottom   = it[slip::n_8c[0]];

  return (center - bottom);
  }
};

template<typename Real>
struct __dz_bck<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dz_bck(_II it)
  {
    Real center = *it;
    Real bottom = it[slip::n_6c[5]];

  return (center - bottom);
  }
};


template<typename Real>
struct __dz_bck<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dz_bck(_II it)
                {
                                Real center = *it;
                                Real bottom = it[slip::n_4d_8c[6]];

                                return (center - bottom);
                }
};

template<typename Real> struct DzBck
{
  
  DzBck(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dz_bck<Real,_Category>::dz_bck(it);
  }
};


template<typename Real,typename>
struct __dt_bck
{
                template <typename _II>
                static Real
                dt_bck(_II it)
                {
                                Real prev    = *(it - 1);
                                Real center  = *it;

                                return (center - prev);
                }
};

template<typename Real>
struct __dt_bck<Real,std::random_access_iterator_tag>
{
                template <typename _II>
                static Real
                dt_bck(_II it)
                {
                                return __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
                }
};

template<typename Real>
struct __dt_bck<Real, std::random_access_iterator2d_tag>
{
                template <typename _II>
                static Real
                dt_bck(_II it)
                {
                                Real center   = *it;
                                Real prev   = it[slip::n_8c[5]];

                                return (center - prev);
                }
};

template<typename Real>
struct __dt_bck<Real, std::random_access_iterator3d_tag>
{
                template <typename _II>
                static Real
                dt_bck(_II it)
                {
                                Real center = *it;
                                Real prev = it[slip::n_26c[5]];

                                return (center - prev);
                }
};

template<typename Real>
struct __dt_bck<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dt_bck(_II it)
                {
                                Real center = *it;
                                Real prev = it[slip::n_4d_8c[0]];

                                return (center - prev);
                }
};

template<typename Real> struct DtBck
{

                DtBck(){}
                template<typename Iterator>
                Real operator() (Iterator& it) const
                {
                                typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
                                return __dt_bck<Real,_Category>::dt_bck(it);
                }
};


//
// Prewittx
//
template<typename Real,typename>
struct __prewittx
{
  template <typename _II>
  static Real
  prewittx(_II it)
  {
    return  __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __prewittx<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  prewittx(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __prewittx<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  prewittx(_II it)
  {
    
    Real NE = it[slip::n_8c[1]];
    Real  E = it[slip::n_8c[0]];
    Real NW = it[slip::n_8c[3]];
    Real  W = it[slip::n_8c[4]];
    Real SE = it[slip::n_8c[7]];
    Real SW = it[slip::n_8c[5]];
    
    return Real(0.5) * (((NE - NW) + (E - W))  + (SE - SW));
  }
};

template<typename Real>
struct __prewittx<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  prewittx(_II it)
  {
     //right
    Real RE  = it[slip::n_26c[1]];
    Real RNE = it[slip::n_26c[2]];
    Real RN  = it[slip::n_26c[10]];
    Real RNW = it[slip::n_26c[19]];
    Real RW  = it[slip::n_26c[18]];
    Real RSW = it[slip::n_26c[17]];
    Real RS  = it[slip::n_26c[25]];
    Real RSE = it[slip::n_26c[8]];
    Real S_RC  = it[slip::n_26c[9]];
   
    //left
    Real LE  = it[slip::n_26c[5]];
    Real LNE = it[slip::n_26c[4]];
    Real LN  = it[slip::n_26c[12]];
    Real LNW = it[slip::n_26c[21]];
    Real LW  = it[slip::n_26c[14]];
    Real LSW = it[slip::n_26c[15]];
    Real LS  = it[slip::n_26c[23]];
    Real LSE = it[slip::n_26c[6]];
    Real LC  = it[slip::n_26c[22]];
    
    return  Real(0.5) * ((RE - LE)  + (RNE - LNE) + (RN - LN)
                                      + (RNW - LNW)+ (RW - LW)   + (RSW - LSW) 
                                      + (RS - LS)  + (RSE - LSE) + (S_RC - LC));
  }
};

template<typename Real> struct Prewittx
{
  
  Prewittx(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __prewittx<Real,_Category>::prewittx(it);
  }
};

//
// Prewitty
//
template<typename Real,typename>
struct __prewitty
{
  template <typename _II>
  static Real
  prewitty(_II it)
  {
    return  __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __prewitty<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  prewitty(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __prewitty<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  prewitty(_II it)
  {
    
    Real NE = it[slip::n_8c[1]];
    Real  N = it[slip::n_8c[2]];
    Real NW = it[slip::n_8c[3]];
    Real  S = it[slip::n_8c[6]];
    Real SE = it[slip::n_8c[7]];
    Real SW = it[slip::n_8c[5]];
    
    return Real(0.5) * (((NE - SE) + (N - S))  + (NW - SW));
  }
};

template<typename Real>
struct __prewitty<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  prewitty(_II it)
  {
     //top
    Real TE  = it[slip::n_26c[10]];
    Real TNE = it[slip::n_26c[2]];
    Real TN  = it[slip::n_26c[3]];
    Real TNW = it[slip::n_26c[4]];
    Real TW  = it[slip::n_26c[12]];
    Real TSW = it[slip::n_26c[21]];
    Real TS  = it[slip::n_26c[20]];
    Real TSE = it[slip::n_26c[19]];
    Real TC  = it[slip::n_26c[11]];
   
    //bottom
    Real BE  = it[slip::n_26c[25]];
    Real BNE = it[slip::n_26c[8]];
    Real BN  = it[slip::n_26c[7]];
    Real BNW = it[slip::n_26c[6]];
    Real BW  = it[slip::n_26c[23]];
    Real BSW = it[slip::n_26c[15]];
    Real BS  = it[slip::n_26c[16]];
    Real BSE = it[slip::n_26c[17]];
    Real S_BC  = it[slip::n_26c[24]];
    
    return Real(0.5) * ((TE - BE)  + (TNE - BNE) + (TN - BN)
                                      + (TNW - BNW)+ (TW - BW)   + (TSW - BSW) 
                                      + (TS - BS)  + (TSE - BSE) + (TC - S_BC));
  }
};

template<typename Real> struct Prewitty
{
  
  Prewitty(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __prewitty<Real,_Category>::prewitty(it);
  }
};

//
// Prewittz
//
template<typename Real,typename>
struct __prewittz
{
  template <typename _II>
  static Real
  prewittz(_II it)
  {
    return  __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __prewittz<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  prewittz(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};


template<typename Real>
struct __prewittz<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  prewittz(_II it)
  {
     //front
    Real FE  = it[slip::n_26c[5]];
    Real FNE = it[slip::n_26c[4]];
    Real FN  = it[slip::n_26c[3]];
    Real FNW = it[slip::n_26c[2]];
    Real FW  = it[slip::n_26c[1]];
    Real FSW = it[slip::n_26c[8]];
    Real FS  = it[slip::n_26c[7]];
    Real FSE = it[slip::n_26c[6]];
    Real FC  = it[slip::n_26c[0]];
   
    //back
    Real BE  = it[slip::n_26c[14]];
    Real BNE = it[slip::n_26c[21]];
    Real BN  = it[slip::n_26c[20]];
    Real BNW = it[slip::n_26c[19]];
    Real BW  = it[slip::n_26c[18]];
    Real BSW = it[slip::n_26c[17]];
    Real BS  = it[slip::n_26c[16]];
    Real BSE = it[slip::n_26c[15]];
    Real S_BC  = it[slip::n_26c[13]];
    
    return Real(0.5) * ((FE - BE)  + (FNE - BNE) + (FN - BN)
                                      + (FNW - BNW)+ (FW - BW)   + (FSW - BSW) 
                                      + (FS - BS)  + (FSE - BSE) + (FC - S_BC));
  }
};

template<typename Real> struct Prewittz
{
  
  Prewittz(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __prewittz<Real,_Category>::prewittz(it);
  }
};

//
// Sobelx
//
template<typename Real,typename>
struct __sobelx
{
  template <typename _II>
  static Real
  sobelx(_II it)
  {
    return  __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __sobelx<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  sobelx(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __sobelx<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  sobelx(_II it)
  {
    
    Real NE = it[slip::n_8c[1]];
    Real  E = it[slip::n_8c[0]];
    Real NW = it[slip::n_8c[3]];
    Real  W = it[slip::n_8c[4]];
    Real SE = it[slip::n_8c[7]];
    Real SW = it[slip::n_8c[5]];
    
    return Real(0.5) * (((NE - NW) + Real(2.0) * (E - W))  + (SE - SW));
  }
};

template<typename Real>
struct __sobelx<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  sobelx(_II it)
  {
    //right
    Real RE  = it[slip::n_26c[1]];
    Real RNE = it[slip::n_26c[2]];
    Real RN  = it[slip::n_26c[10]];
    Real RNW = it[slip::n_26c[19]];
    Real RW  = it[slip::n_26c[18]];
    Real RSW = it[slip::n_26c[17]];
    Real RS  = it[slip::n_26c[25]];
    Real RSE = it[slip::n_26c[8]];
    Real S_RC  = it[slip::n_26c[9]];
   
    //left
    Real LE  = it[slip::n_26c[5]];
    Real LNE = it[slip::n_26c[4]];
    Real LN  = it[slip::n_26c[12]];
    Real LNW = it[slip::n_26c[21]];
    Real LW  = it[slip::n_26c[14]];
    Real LSW = it[slip::n_26c[15]];
    Real LS  = it[slip::n_26c[23]];
    Real LSE = it[slip::n_26c[6]];
    Real LC  = it[slip::n_26c[22]];
    
    return (Real(2.0) * (RE - LE) + (RNE - LNE) + Real(2.0) * (RN - LN)
                    + (RNW - LNW) + Real(2.0) * (RW - LW) + (RSW - LSW) + Real(2.0) * (RS - LS) + (RSE - LSE)+ Real(4.0) * (S_RC - LC));
  }
};

template<typename Real> struct Sobelx
{
  
  Sobelx(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __sobelx<Real,_Category>::sobelx(it);
  }
};


//
// Sobely
//
template<typename Real,typename>
struct __sobely
{
  template <typename _II>
  static Real
  sobely(_II it)
  {
    return  __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __sobely<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  sobely(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __sobely<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  sobely(_II it)
  {
    
    Real NE = it[slip::n_8c[1]];
    Real  N = it[slip::n_8c[2]];
    Real NW = it[slip::n_8c[3]];
    Real  S = it[slip::n_8c[6]];
    Real SE = it[slip::n_8c[7]];
    Real SW = it[slip::n_8c[5]];
    
    return Real(0.5) * (((NE - SE) + Real(2.0) * (N - S))  + (NW - SW));
  }
};

template<typename Real>
struct __sobely<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  sobely(_II it)
  {
      //top
    Real TE  = it[slip::n_26c[10]];
    Real TNE = it[slip::n_26c[2]];
    Real TN  = it[slip::n_26c[3]];
    Real TNW = it[slip::n_26c[4]];
    Real TW  = it[slip::n_26c[12]];
    Real TSW = it[slip::n_26c[21]];
    Real TS  = it[slip::n_26c[20]];
    Real TSE = it[slip::n_26c[19]];
    Real TC  = it[slip::n_26c[11]];
   
    //bottom
    Real BE  = it[slip::n_26c[25]];
    Real BNE = it[slip::n_26c[8]];
    Real BN  = it[slip::n_26c[7]];
    Real BNW = it[slip::n_26c[6]];
    Real BW  = it[slip::n_26c[23]];
    Real BSW = it[slip::n_26c[15]];
    Real BS  = it[slip::n_26c[16]];
    Real BSE = it[slip::n_26c[17]];
    Real S_BC  = it[slip::n_26c[24]];
    
    return (Real(2.0) * (TE - BE) + (TNE - BNE) + Real(2.0) * (TN - BN)
                    + (TNW - BNW) + Real(2.0) * (TW - BW) + (TSW - BSW) + Real(2.0) * (TS - BS) + (TSE - BSE)+ Real(4.0) * (TC - S_BC));
  }
};

template<typename Real> struct Sobely
{
  
  Sobely(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __sobely<Real,_Category>::sobely(it);
  }
};

//
// Sobelz
//
template<typename Real,typename>
struct __sobelz
{
  template <typename _II>
  static Real
  sobelz(_II it)
  {
    return  __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __sobelz<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  sobelz(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};



template<typename Real>
struct __sobelz<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  sobelz(_II it)
  {
   
      //front
    Real FE  = it[slip::n_26c[5]];
    Real FNE = it[slip::n_26c[4]];
    Real FN  = it[slip::n_26c[3]];
    Real FNW = it[slip::n_26c[2]];
    Real FW  = it[slip::n_26c[1]];
    Real FSW = it[slip::n_26c[8]];
    Real FS  = it[slip::n_26c[7]];
    Real FSE = it[slip::n_26c[6]];
    Real FC  = it[slip::n_26c[0]];
   
    //back
    Real BE  = it[slip::n_26c[14]];
    Real BNE = it[slip::n_26c[21]];
    Real BN  = it[slip::n_26c[20]];
    Real BNW = it[slip::n_26c[19]];
    Real BW  = it[slip::n_26c[18]];
    Real BSW = it[slip::n_26c[17]];
    Real BS  = it[slip::n_26c[16]];
    Real BSE = it[slip::n_26c[15]];
    Real S_BC  = it[slip::n_26c[13]];

    return (Real(2.0) * (FE - BE) + (FNE - BNE) + Real(2.0) * (FN - BN)
                    + (FNW - BNW) + Real(2.0) * (FW - BW) + (FSW - BSW) + Real(2.0) * (FS - BS) + (FSE - BSE)+ Real(4.0) * (FC - S_BC));
  }
};

template<typename Real> struct Sobelz
{
  
  Sobelz(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __sobelz<Real,_Category>::sobelz(it);
  }
};

//
//dv first derivative in the direction of v
//
template<typename Real,typename>
struct __dv
{
  template <typename _II, typename Block>
  static Real
  dv(_II it, Block v)
  {
    assert(Block::SIZE == 1);
    Real dx = __dx<Real,std::random_access_iterator_tag>::dx(it);
    if(v[0] < Real(0))
      {
                dx = -dx;
      }
    return dx;
  }
};

template<typename Real>
struct __dv<Real,std::random_access_iterator_tag>
{
  template <typename _II, typename Block>
   static Real
  dv(_II it, Block v)
  {
    assert(Block::SIZE == 1);
    return __dv<Real,std::random_access_iterator_tag>::dv(it,v);
  }
};

template<typename Real>
struct __dv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II, typename Block>
   static Real
  dv(_II it, Block v)
  {
    assert(Block::SIZE == 2);
    Real result = Real(0);
    Real norm = slip::Euclidean_norm<Real>(v.begin(),v.end());

    if(norm != Real(0))
      {
                Block grad;
                grad[0] = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
                grad[1] = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
                result = std::inner_product(grad.begin(),grad.end(),v.begin(),Real()) / norm;
      }
    return result;
  }
};

template<typename Real>
struct __dv<Real, std::random_access_iterator3d_tag>
{
  template <typename _II, typename Block>
   static Real
  dv(_II it, Block v)
  {
    assert(Block::SIZE == 3);
    Real result = Real(0);
    Real norm = slip::Euclidean_norm<Real>(v.begin(),v.end());
    if(norm != Real(0))
      {
                Block grad;
                grad[0] = __dx<Real,std::random_access_iterator3d_tag>::dx(it);
                grad[1] = __dy<Real,std::random_access_iterator3d_tag>::dy(it);
                grad[2] = __dz<Real,std::random_access_iterator3d_tag>::dz(it);
                result = std::inner_product(grad.begin(),grad.end(),v.begin(),Real()) / norm;
      }
    return result;
  }
};


template<typename Real>
struct __dv<Real, std::random_access_iterator4d_tag>
{
                template <typename _II, typename Block>
                static Real
                dv(_II it, Block v)
                {
                                assert(Block::SIZE == 4);
                                Real result = Real(0);
                                Real norm = slip::Euclidean_norm<Real>(v.begin(),v.end());
                                if(norm != Real(0))
                                {
                                                Block grad;
                                                grad[0] = __dt<Real,std::random_access_iterator4d_tag>::dt(it);
                                                grad[1] = __dx<Real,std::random_access_iterator4d_tag>::dx(it);
                                                grad[2] = __dy<Real,std::random_access_iterator4d_tag>::dy(it);
                                                grad[3] = __dz<Real,std::random_access_iterator4d_tag>::dz(it);
                                                result = std::inner_product(grad.begin(),grad.end(),v.begin(),Real()) / norm;
                                }
                                return result;
                }
};


template<typename Real, typename Block> struct Dv 
{
  Dv(Block v)
    :v_(v)
  {}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dv<Real,_Category>::dv(it,v_);
  }
  Block v_;
};


//
// gradient norm
//
template<typename Real,typename>
struct __grad
{
  template <typename _II>
  static Real
  grad(_II it)
  {
    Real dx = __dx<Real,std::random_access_iterator_tag>::dx(it);
    
    return std::abs(dx);
  }
};

template<typename Real>
struct __grad<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  grad(_II it)
  {
    Real dx = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real dy = __dy<Real,std::random_access_iterator2d_tag>::dy(it);

    return std::sqrt((dx * dx) + (dy * dy));
    
  }
};

template<typename Real>
struct __grad<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  grad(_II it)
  {
    Real dx = __dx<Real,std::random_access_iterator3d_tag>::dx(it);
    Real dy = __dy<Real,std::random_access_iterator3d_tag>::dy(it);
    Real dz = __dz<Real,std::random_access_iterator3d_tag>::dz(it);

    return std::sqrt((dx * dx) + (dy * dy) + (dz * dz));
    
  }
};



template<typename Real>
struct __grad<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                grad(_II it)
                {
                                Real dt = __dt<Real,std::random_access_iterator4d_tag>::dt(it);
                                Real dx = __dx<Real,std::random_access_iterator4d_tag>::dx(it);
                                Real dy = __dy<Real,std::random_access_iterator4d_tag>::dy(it);
                                Real dz = __dz<Real,std::random_access_iterator4d_tag>::dz(it);

                                return std::sqrt((dt * dt) + (dx * dx) + (dy * dy) + (dz * dz));

                }
};


template<typename Real> struct Grad
{ 
  Grad(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __grad<Real,_Category>::grad(it);
  }
};




//
// dxx local derivatives
//

template<typename Real,typename>
struct __dxx
{
  template <typename _II>
  static Real
  dxx(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxx<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dxx(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxx<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dxx(_II it)
  {
    Real left    = it[slip::n_8c[4]];
    Real center  = *it;
    Real right   = it[slip::n_8c[0]];

    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxx<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dxx(_II it)
  {
    Real left    = it[slip::n_6c[3]];
    Real center  = *it;
    Real right   = it[slip::n_6c[1]];

    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxx<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dxx(_II it)
                {
                                Real left    = it[slip::n_4d_8c[4]];
                                Real center  = *it;
                                Real right   = it[slip::n_4d_8c[2]];

                                return ((right + left) - Real(2) * center);
                }
};


template<typename Real> struct Dxx
{
 
  Dxx(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dxx<Real,_Category>::dxx(it);
  }
};


//
// dyy local derivatives
//
template<typename Real,typename>
struct __dyy
{
  template <typename _II>
  static Real
  dyy(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dyy<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dyy(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __dyy<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dyy(_II it)
  {
    Real top      = it[slip::n_8c[2]];
    Real center   = *it;
    Real bottom   = it[slip::n_8c[6]];

  return ((top + bottom) - Real(2) * center);
  }
};

template<typename Real>
struct __dyy<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dyy(_II it)
  {
    Real top    = it[slip::n_6c[2]];
    Real center   = *it;
    Real bottom = it[slip::n_6c[4]];

  return ((top + bottom) - Real(2) * center);
  }
};

template<typename Real>
struct __dyy<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dyy(_II it)
                {
                                Real top    = it[slip::n_4d_8c[3]];
                                Real center   = *it;
                                Real bottom = it[slip::n_4d_8c[5]];

                                return ((top + bottom) - Real(2) * center);
                }
};

template<typename Real> struct Dyy
{
 
  Dyy(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dyy<Real,_Category>::dyy(it);
  }
};

//
// dzz local derivatives
//
template<typename Real,typename>
struct __dzz
{
  template <typename _II>
  static Real
  dzz(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dzz<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dzz(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __dzz<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dzz(_II it)
  {
    Real top      = it[slip::n_8c[5]];
    Real center   = *it;
    Real bottom   = it[slip::n_8c[0]];

  return ((top + bottom) - Real(2) * center);
  }
};

template<typename Real>
struct __dzz<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dzz(_II it)
  {
    Real top    = it[slip::n_6c[5]];
    Real center = *it;
    Real bottom = it[slip::n_6c[0]];

  return ((top + bottom) - Real(2) * center);
  }
};


template<typename Real>
struct __dzz<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dzz(_II it)
                {
                                Real top    = it[slip::n_4d_8c[6]];
                                Real center = *it;
                                Real bottom = it[slip::n_4d_8c[1]];

                                return ((top + bottom) - Real(2) * center);
                }
};

template<typename Real> struct Dzz
{
 
  Dzz(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dzz<Real,_Category>::dzz(it);
  }
};

template<typename Real,typename>
struct __dtt
{
                template <typename _II>
                static Real
                dtt(_II it)
                {
                                Real prev   = *(it - 1);
                                Real center  = *it;
                                Real next   = *(it + 1);

                                return ((next + prev) - Real(2) * center);
                }
};

template<typename Real>
struct __dtt<Real,std::random_access_iterator_tag>
{
                template <typename _II>
                static Real
                dtt(_II it)
                {
                                return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
                }
};

template<typename Real>
struct __dtt<Real, std::random_access_iterator2d_tag>
{
                template <typename _II>
                static Real
                dtt(_II it)
                {
                                Real prev     = it[slip::n_8c[5]];
                                Real center   = *it;
                                Real next   = it[slip::n_8c[0]];

                                return ((next + prev) - Real(2) * center);
                }
};

template<typename Real>
struct __dtt<Real, std::random_access_iterator3d_tag>
{
                template <typename _II>
                static Real
                dtt(_II it)
                {
                                Real prev    = it[slip::n_26c[5]];
                                Real center = *it;
                                Real next = it[slip::n_26c[19]];

                                return ((next + prev) - Real(2) * center);
                }
};

template<typename Real>
struct __dtt<Real, std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                dtt(_II it)
                {
                                Real prev    = it[slip::n_4d_8c[0]];
                                Real center = *it;
                                Real next = it[slip::n_4d_8c[7]];

                                return ((next + prev) - Real(2) * center);
                }
};

template<typename Real> struct Dtt
{

                Dtt(){}
                template<typename Iterator>
                Real operator() (Iterator& it) const
                {
                                typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
                                return __dtt<Real,_Category>::dtt(it);
                }
};

template<typename Real,typename>
struct __dxy
{
  template <typename _II>
  static Real
  dxy(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxy<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dxy(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxy<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dxy(_II it)
  {
    
    Real NE = it[slip::n_8c[1]];
    Real NW = it[slip::n_8c[3]];
    Real SE = it[slip::n_8c[7]];
    Real SW = it[slip::n_8c[5]];
    
    return ((SW + NE) - (SE + NW)) * Real(0.25);
  }
};

template<typename Real>
struct __dxy<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dxy(_II it)
  {
    Real NE = it[slip::n_26c[10]];
    Real NW = it[slip::n_26c[12]];
    Real SE = it[slip::n_26c[25]];
    Real SW = it[slip::n_26c[23]];
  
    return ((SW + NE) - (SE + NW)) * Real(0.25);
  }
};

template<typename Real> struct Dxy
{
  
  Dxy(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dxy<Real,_Category>::dxy(it);
  }
};

//
//dxz
template<typename Real,typename>
struct __dxz
{
  template <typename _II>
  static Real
  dxz(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __dxz<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dxz(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};


template<typename Real>
struct __dxz<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dxz(_II it)
  {
    Real NE = it[slip::n_26c[1]];
    Real NW = it[slip::n_26c[5]];
    Real SE = it[slip::n_26c[18]];
    Real SW = it[slip::n_26c[14]];
  
    return ((SW + NE) - (SE + NW)) * Real(0.25);
  }
};

template<typename Real> struct Dxz
{
  
  Dxz(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dxz<Real,_Category>::dxz(it);
  }
};

//
//dyz
template<typename Real,typename>
struct __dyz
{
  template <typename _II>
  static Real
  dyz(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __dyz<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dyz(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};



template<typename Real>
struct __dyz<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dyz(_II it)
  {
    Real NE = it[slip::n_26c[3]];
    Real NW = it[slip::n_26c[20]];
    Real SE = it[slip::n_26c[18]];
    Real SW = it[slip::n_26c[7]];
  
    return ((SW + NE) - (SE + NW)) * Real(0.25);
  }
};

template<typename Real> struct Dyz
{
  
  Dyz(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dyz<Real,_Category>::dyz(it);
  }
};



//
//dxy Weickert
//

template<typename Real,typename>
struct __dxy_w
{
  template <typename _II>
  static Real
  dxy_w(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxy_w<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  dxy_w(_II it)
  {
    Real left    = *(it - 1);
    Real center  = *it;
    Real right   = *(it + 1);
    
    return ((right + left) - Real(2) * center);
  }
};

template<typename Real>
struct __dxy_w<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  dxy_w(_II it)
  {
    Real gx  = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy  = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gxgy = gx * gy;
    
    Real dxy = Real(0);
    
    Real center = *it;
    Real E  = it[slip::n_8c[0]];
    Real NE = it[slip::n_8c[1]];
    Real N  = it[slip::n_8c[2]];
    Real NW = it[slip::n_8c[3]];
    Real W  = it[slip::n_8c[4]];
    Real SW = it[slip::n_8c[5]];
    Real S  = it[slip::n_8c[6]];
    Real SE = it[slip::n_8c[7]];
    if(gxgy > 0)
      {
                dxy = (SW + NE + Real(2) * center - E - S - N - W) * Real(0.5);
      }
    else
      {
                dxy = (E + S + N + W - NW - SE - Real(2) * center) * Real(0.5);
      }

    return dxy;
  }
};

template<typename Real>
struct __dxy_w<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  dxy_w(_II it)
  {
    Real gx  = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy  = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gxgy = gx * gy;
    
    Real dxy = Real(0);
    
    Real center = *it;
    Real E  = it[slip::n_26c[9]];
    Real NE = it[slip::n_26c[10]];
    Real N  = it[slip::n_26c[11]];
    Real NW = it[slip::n_26c[12]];
    Real W  = it[slip::n_26c[22]];
    Real SW = it[slip::n_26c[23]];
    Real S  = it[slip::n_26c[24]];
    Real SE = it[slip::n_26c[25]];
 
    if(gxgy > 0)
      {
                dxy = (SW + NE + Real(2) * center - E - S - N - W) * Real(0.5);
      }
    else
      {
                dxy = (E + S + N + W - NW - SE - Real(2) * center) * Real(0.5);
      }

    return dxy;
   
  }
};

template<typename Real> struct Dxy_W
{
  
  Dxy_W(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dxy_w<Real,_Category>::dxy_w(it);
  }
};

//
// divergence local derivatives
//
template<typename Real,typename>
struct __div
{
  template <typename _II>
  static Real
  div(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __div<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  div(_II it)
  {
    return __dx<Real,std::random_access_iterator_tag>::dx(it);
  }
};

template<typename Real>
struct __div<Real,std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  div(_II it)
  {
    return __dx<Real,std::random_access_iterator2d_tag>::dx(it)
          +__dy<Real,std::random_access_iterator2d_tag>::dy(it);
  }
};

template<typename Real>
struct __div<Real,std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  div(_II it)
  {
    return __dx<Real,std::random_access_iterator3d_tag>::dx(it)
          +__dy<Real,std::random_access_iterator3d_tag>::dy(it)
          +__dz<Real,std::random_access_iterator3d_tag>::dz(it);
  }
};

template<typename Real>
struct __div<Real,std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                div(_II it)
                {
                                return __dt<Real,std::random_access_iterator4d_tag>::dt(it)
                                                                +__dx<Real,std::random_access_iterator4d_tag>::dx(it)
                                                                +__dy<Real,std::random_access_iterator4d_tag>::dy(it)
                                                                +__dz<Real,std::random_access_iterator4d_tag>::dz(it);
                }
};


template<typename Real> struct Div
{
  
  Div(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __div<Real,_Category>::div(it);
  }
};




//
//  local laplacian
//
template<typename Real,typename>
struct __lap
{
  template <typename _II>
  static Real
  lap(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __lap<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  lap(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __lap<Real,std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  lap(_II it)
  {
    return __dxx<Real,std::random_access_iterator2d_tag>::dxx(it)
          +__dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
  }
};

template<typename Real>
struct __lap<Real,std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  lap(_II it)
  {
    return __dxx<Real,std::random_access_iterator3d_tag>::dxx(it)
          +__dyy<Real,std::random_access_iterator3d_tag>::dyy(it)
          +__dzz<Real,std::random_access_iterator3d_tag>::dzz(it);
  }
};

template<typename Real>
struct __lap<Real,std::random_access_iterator4d_tag>
{
                template <typename _II>
                static Real
                lap(_II it)
                {
                                return __dtt<Real,std::random_access_iterator4d_tag>::dtt(it)
                                                                +__dxx<Real,std::random_access_iterator4d_tag>::dxx(it)
                                                                +__dyy<Real,std::random_access_iterator4d_tag>::dyy(it)
                                                                +__dzz<Real,std::random_access_iterator4d_tag>::dzz(it);
                }
};


template<typename Real> struct Lap 
{
  
  Lap(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __lap<Real,_Category>::lap(it);
  }
};

template<typename Real,typename>
struct __lap_8c
{
  template <typename _II>
  static Real
  lap(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};
template<typename Real>
struct __lap_8c<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  lap(_II it)
  {
    
    const Real gamma = Real(2) / Real(3);
    Real gamma2 = Real(1) - gamma;
    Real inv_2h2 = Real(0.5);

    Real Ixx  = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real Iyy  = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real NE = it[slip::n_8c[1]];
    Real NW = it[slip::n_8c[3]];
    Real SE = it[slip::n_8c[7]];
    Real SW = it[slip::n_8c[5]];
    Real diag = (((SE + NW) + (SW + NE)) - ((Real(4) * *it))) * inv_2h2;

    return ((gamma2 * Ixx) + (gamma2 * Iyy)) + (gamma * diag);
  }
};
template<typename Real> struct Lap8
{
 
  Lap8(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __lap_8c<Real,_Category>::lap(it);
  }
};


template<typename Real> struct LapLindeberg
{
 
  LapLindeberg(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return slip::constants<Real>::half()*__lap_8c<Real,_Category>::lap(it);
  }
};
//
//  local second derivative in the direction of the gradient
//
template<typename Real,typename>
struct __iee_sapiro
{
  template <typename _II>
  static Real
  iee_sapiro(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __iee_sapiro<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  iee_sapiro(_II it)
  {
    
    Real gx  = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy  = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2 = gx * gx;
    Real gy2 = gy * gy;
    Real g2  = gx2 + gy2;
  
    Real iee = Real(0);

    if(g2 != 0)
      {
                Real lambda1 = gx2;
                Real lambda2 = gy2;
                Real lambda3 = Real(0.5) * (gx * gy);
                Real lambda4 = - lambda3;
                Real lambda0 = Real(-2.0) * (lambda1 + lambda2);
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];

                
                iee = ( lambda4   * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (S + N)
                                + lambda1 * (W + E)
                                + lambda0 * center) / g2 ;
    }
 

    return(iee);

  }
};

template<typename Real>
struct __iee_sapiro<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
  static Real
  iee_sapiro(_II it)
  {
    
    Real dx  = __dx<Real,std::random_access_iterator3d_tag>::dx(it);
    Real dy  = __dy<Real,std::random_access_iterator3d_tag>::dy(it);
    Real dz  = __dz<Real,std::random_access_iterator3d_tag>::dz(it);
    
    Real dx2 = dx * dx;
    Real dy2 = dy * dy;
    Real dz2 = dz + dz;

    Real g2 = (dx2 + dy2) + dz2;

    Real iee = Real(0);

    if(g2 != 0)
      {

                Real dxx = __dxx<Real,std::random_access_iterator3d_tag>::dxx(it);
                Real dyy = __dyy<Real,std::random_access_iterator3d_tag>::dyy(it);
                Real dzz = __dzz<Real,std::random_access_iterator3d_tag>::dzz(it);
                
                Real dxy = __dxy<Real,std::random_access_iterator3d_tag>::dxy(it);
                Real dxz = __dxz<Real,std::random_access_iterator3d_tag>::dxz(it);
                Real dyz = __dyz<Real,std::random_access_iterator3d_tag>::dyz(it);
                
                Real dxdy = dx * dy;
                Real dxdz = dx * dz;
                Real dydz = dy * dz;
                
                Real a = (dy2 + dz2) * dxx;
                Real b = (dx2 + dz2) * dyy;
                Real c = (dx2 + dy2) * dzz;
                Real d = Real(2.0) * (dxdy * dxy);
                Real e = Real(2.0) * (dxdz * dxz);
                Real f = Real(2.0) * (dydz * dyz);

                iee = ((a + b) + (c + d) + (e + f)) / g2;
    }

    return(iee);

  }
};

template<typename Real> struct Iee_Sapiro
{
 
  Iee_Sapiro(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __iee_sapiro<Real,_Category>::iee_sapiro(it);
  }
};

template<typename Real,typename>
struct __iee_lucido
{
  template <typename _II>
  static Real
  iee_lucido(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __iee_lucido<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  iee_lucido(_II it)
  {
    
    Real gx   = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy   = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2  = gx * gx;
    Real gy2  = gy * gy;
    Real g2   = gx2 + gy2;
    Real gxgy = gx * gy;         
    Real iee  = Real(0);

    if(g2 != 0)
      { 
                if(gxgy >= 0)
                  {
                  
      
                Real lambda1 = gx2 - gxgy;
                Real lambda2 = gy2 - gxgy;
                Real lambda3 = gxgy;
                //              Real lambda4 = - lambda3;
                Real lambda0 = Real(-2.0) * ((lambda1 + lambda2) + lambda3);
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];

                iee = ( lambda3 * (NE + SW)
                      + lambda2 * (S + N)
                      + lambda1 * (W + E)
                      + lambda0 * center)/ g2;
                  }
                else
                  {
                    Real lambda1 = gx2  + gxgy;
                    Real lambda2 = gy2  + gxgy;
                    Real lambda4 = - gxgy;
                    Real lambda0 = Real(-2.0) * ((lambda1 + lambda2) + lambda4);
                    Real center = *it;
                    Real E  = it[slip::n_8c[0]];
                    Real N  = it[slip::n_8c[2]];
                    Real NW = it[slip::n_8c[3]];
                    Real W  = it[slip::n_8c[4]];
                    Real S  = it[slip::n_8c[6]];
                    Real SE = it[slip::n_8c[7]];
                    iee = ( lambda4 * (NW + SE)
                                  + lambda2 * (S + N)
                                  + lambda1 * (W + E)
                                  + lambda0 * center)/ g2;
                  }

      }
 

    return(iee);

  }
};
template<typename Real> struct Iee_Lucido
{
  
  Iee_Lucido(){}
  template<typename Iterator>
  Real operator() (Iterator& it)
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __iee_lucido<Real,_Category>::iee_lucido(it);
  }
};
template<typename Real,typename>
struct __iee_augereau
{
  template <typename _II>
  static Real
  iee_augereau(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __iee_augereau<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  iee_augereau(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
    Real absgx = std::abs(gx);
    Real absgy = std::abs(gy);
    Real iee   = Real(0);

    if(g2 != 0)
      { 
                if(absgx >= absgy)
                  {
                    Real lambda0 = Real(-2.0) * gx2;
                    Real lambda1 = gx2 - gy2;
                    Real lambda3 = Real(0.5) * (gy2 + gxgy);
                    Real lambda4 = Real(0.5) * (gy2 - gxgy);
                    
                    Real center = *it;
                    Real E  = it[slip::n_8c[0]];
                    Real NE = it[slip::n_8c[1]];
                    Real NW = it[slip::n_8c[3]];
                    Real W  = it[slip::n_8c[4]];
                    Real SW = it[slip::n_8c[5]];
                    Real SE = it[slip::n_8c[7]];
                
                    iee = ( lambda4 * (NW + SE)
                                    + lambda3 * (NE + SW)
                                    + lambda1 * (W + E)
                                    + lambda0 * center)/ g2;
                  }
                else
                  {
                    Real lambda0 = Real(-2.0) * gy2;
                    Real lambda2 = gy2 - gx2;
                    Real lambda3 = Real(0.5) * (gx2 + gxgy);
                    Real lambda4 = Real(0.5) * (gx2 - gxgy);
                   
                    Real center = *it;
                    Real NE = it[slip::n_8c[1]];
                    Real N  = it[slip::n_8c[2]];
                    Real NW = it[slip::n_8c[3]];
                    Real SW = it[slip::n_8c[5]];
                    Real S  = it[slip::n_8c[6]];
                    Real SE = it[slip::n_8c[7]];
                    iee = ( lambda4 * (NW + SE)
                                  + lambda3 * (NE + SW)
                                  + lambda2 * (S + N)
                                  + lambda0 * center)/ g2;
                  }

      }
 

    return(iee);

  }
};
template<typename Real> struct Iee_Augereau 
{
 
  Iee_Augereau(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __iee_augereau<Real,_Category>::iee_augereau(it);
  }
};




//
template<typename Real,typename>
struct __iee_alvarez
{
  template <typename _II>
  static Real
  iee_alvarez(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __iee_alvarez<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  iee_alvarez(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
   
    Real iee   = Real(0);

    if(g2 != 0)
      { 

                Real lambda0 = Real(0.25);
                Real lambda1 = Real(2.0) * lambda0 - gy2 / g2;
                Real lambda2 = Real(2.0) * lambda0 - gx2 / g2;
                Real lambda3 = -lambda0 + Real(0.5) * ( Real(1) + gxgy / g2);
                Real lambda4 = -lambda0 + Real(0.5) * ( Real(1) - gxgy / g2);
                lambda0 = Real(-1.0);//-4.0*0.25 
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];
                
                iee = (   lambda4 * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (N + S)
                                + lambda1 * (W + E)
                                + lambda0 * center);
      }
  
    return(iee);

  }
};
template<typename Real> struct Iee_Alvarez
{
  
  Iee_Alvarez(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __iee_alvarez<Real,_Category>::iee_alvarez(it);
  }
};

//
template<typename Real,typename>
struct __iee_alvarez2
{
  template <typename _II>
  static Real
  iee_alvarez2(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __iee_alvarez2<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  iee_alvarez2(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
   
    Real iee   = Real(0);

    if(g2 != 0)
      { 
                Real absgx = std::abs(gx);
                Real absgy = std::abs(gy);
                Real lambda0 = Real(0);
                if ( (absgx > absgy) && (gxgy>=0) )
                {
                  lambda0 = ( Real(2.0) * gx2 + gy2 - gxgy ) / Real(4.0);
                
                } 
      else if ( (absgx <= absgy) && (gxgy >= 0) )
                {
                  lambda0 = ( Real(2.0) * gy2 + gx2 - gxgy ) / Real(4.0);
                 
                }
      else if ( (absgx > absgy) && (gxgy <= 0) )
                {
                  
                  lambda0 = ( Real(2.0) * gx2 + gy2 + gxgy ) / Real(4.0);       
                 
                }    
      else
                {
                  lambda0 = ( gy2 + Real(2.0) * gx2 - gxgy ) / Real(4.0);
                
                } 

                Real lambda1 = Real(2.0) * lambda0 - gy2;
                Real lambda2  = Real(2.0) *lambda0 - gx2;
                Real lambda3 = -lambda0 + Real(0.5) * ( gx2 + gxgy + gy2 );
                Real lambda4 = -lambda0 + Real(0.5) * ( gx2 - gxgy + gy2 );
                lambda0 = -4.0*lambda0;

                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];
                
                iee = (   lambda4 * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (N + S)
                                + lambda1 * (W + E)
                                + lambda0 * center)/ g2;
      }
  
    return(iee);

  }
};
template<typename Real> struct Iee_Alvarez2 
{
  Iee_Alvarez2(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __iee_alvarez2<Real,_Category>::iee_alvarez2(it);
  }
};
//
template<typename Real,typename>
struct __iee_cohignac
{
  template <typename _II>
  static Real
  iee_cohignac(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __iee_cohignac<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  iee_cohignac(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
   
    Real iee   = Real(0);

    if(g2 != 0)
      { 
                Real lambda0 = Real(0.5) - gx2/g2 + (gx2 / g2) * (gy2 / g2);
                Real lambda1 = Real(2.0) * lambda0 - gy2 /g2;
                Real lambda2 = Real(2.0) * lambda0 - gx2 / g2;
                Real lambda3 = -lambda0 + Real(0.5) * (Real(1.0) + gxgy / g2);
                Real lambda4 = -lambda0 + Real(0.5) * (Real(1.0) - gxgy / g2);
                lambda0 = Real(-4.0) * lambda0;
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];
                
                iee = (   lambda4 * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (N + S)
                                + lambda1 * (W + E)
                                + lambda0 * center);
      }
  
    return(iee);

  }
};
template<typename Real> struct Iee_Cohignac
{
  Iee_Cohignac(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __iee_cohignac<Real,_Category>::iee_cohignac(it);
  }
};


//
//  local second derivative in the direction 
//  perpendicular to the gradient
//
template<typename Real,typename>
struct __inn_sapiro
{
  template <typename _II>
  static Real
  inn_sapiro(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __inn_sapiro<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  inn_sapiro(_II it)
  {
    
    Real gx  = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy  = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2 = gx * gx;
    Real gy2 = gy * gy;
    Real g2  = gx2 + gy2;
    Real gxgy = gx * gy;         
    Real inn = Real(0);

    if(g2 != 0)
      {
                Real lambda1 = gy2;
                Real lambda2 = gx2;
                Real lambda4 = Real(0.5) * gxgy;
                Real lambda3 = - lambda4;
                Real lambda0 = Real(-2.0) * (lambda1 + lambda2);
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];

                
                inn = ( lambda4   * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (S + N)
                                + lambda1 * (W + E)
                                + lambda0 * center) / g2 ;
    }
 

    return(inn);

  }
};

template<typename Real>
struct __inn_sapiro<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  inn_sapiro(_II it)
  {
    
    Real dx  = __dx<Real,std::random_access_iterator3d_tag>::dx(it);
    Real dy  = __dy<Real,std::random_access_iterator3d_tag>::dy(it);
    Real dz  = __dz<Real,std::random_access_iterator3d_tag>::dz(it);
    
    Real dx2 = dx * dx;
    Real dy2 = dy * dy;
    Real dz2 = dz + dz;

    Real g2 = (dx2 + dy2) + dz2;

    Real inn = Real(0);

    if(g2 != 0)
      {

                Real dxx = __dxx<Real,std::random_access_iterator3d_tag>::dxx(it);
                Real dyy = __dyy<Real,std::random_access_iterator3d_tag>::dyy(it);
                Real dzz = __dzz<Real,std::random_access_iterator3d_tag>::dzz(it);
                
                Real dxy = __dxy<Real,std::random_access_iterator3d_tag>::dxy(it);
                Real dxz = __dxz<Real,std::random_access_iterator3d_tag>::dxz(it);
                Real dyz = __dyz<Real,std::random_access_iterator3d_tag>::dyz(it);
                
                Real dxdy = dx * dy;
                Real dxdz = dx * dz;
                Real dydz = dy * dz;
                
                Real a = (dy2 + dz2) * dxx;
                Real b = (dx2 + dz2) * dyy;
                Real c = (dx2 + dy2) * dzz;
                Real d = - Real(2.0) * (dxdy * dxy);
                Real e = - Real(2.0) * (dxdz * dxz);
                Real f = - Real(2.0) * (dydz * dyz);

                inn = ((a + b) + (c + d) + (e + f)) / g2;
    }
     return(inn);

  }
};

template<typename Real> struct Inn_Sapiro
{
  Inn_Sapiro(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __inn_sapiro<Real,_Category>::inn_sapiro(it);
  }
};

template<typename Real,typename>
struct __inn_lucido
{
  template <typename _II>
  static Real
  inn_lucido(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __inn_lucido<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  inn_lucido(_II it)
  {
    
    Real gx   = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy   = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2  = gx * gx;
    Real gy2  = gy * gy;
    Real g2   = gx2 + gy2;
    Real gxgy = gx * gy;         
    Real inn  = Real(0);

    if(g2 != 0)
      { 
                if(gxgy >= 0)
                  {
                  
      
                Real lambda1 = gy2 - gxgy;
                Real lambda2 = gx2 - gxgy;
                Real lambda4 = gxgy;
                Real lambda0 = Real(-2.0) * ((lambda1 + lambda2) + lambda4);
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];

                inn = ( lambda4 * (NW + SE)
                      + lambda2 * (S + N)
                      + lambda1 * (W + E)
                      + lambda0 * center)/ g2;
                  }
                else
                  {
                    Real lambda1 = gy2  + gxgy;
                    Real lambda2 = gx2  + gxgy;
                    Real lambda3 = - gxgy;
                    Real lambda0 = Real(-2.0) * ((lambda1 + lambda2) + lambda3);
                    Real center = *it;
                    Real E  = it[slip::n_8c[0]];
                    Real N  = it[slip::n_8c[2]];
                    Real SW = it[slip::n_8c[5]];
                    Real W  = it[slip::n_8c[4]];
                    Real S  = it[slip::n_8c[6]];
                    Real NE = it[slip::n_8c[1]];
                    inn = ( lambda3 * (SW + NE)
                                  + lambda2 * (S + N)
                                  + lambda1 * (W + E)
                                  + lambda0 * center)/ g2;
                  }

      }
 

    return(inn);

  }
};

template<typename Real> struct Inn_Lucido 
{
  Inn_Lucido(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __inn_lucido<Real,_Category>::inn_lucido(it);
  }
};

//
template<typename Real,typename>
struct __inn_augereau
{
  template <typename _II>
  static Real
  inn_augereau(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __inn_augereau<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  inn_augereau(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
    Real absgx = std::abs(gx);
    Real absgy = std::abs(gy);
    Real inn   = Real(0);

    if(g2 != 0)
      { 
                if(absgx >= absgy)
                  {
                    Real lambda0 = Real(-2.0) * gx2;
                    Real lambda2 = gx2 - gy2;
                    Real lambda3 = Real(0.5) * (gy2 - gxgy);
                    Real lambda4 = Real(0.5) * (gy2 + gxgy);
                    
                    Real center = *it;
                   
                    Real NE = it[slip::n_8c[1]];
                    Real N  = it[slip::n_8c[2]];
                    Real NW = it[slip::n_8c[3]];
                    Real SW = it[slip::n_8c[5]];
                    Real S  = it[slip::n_8c[6]];
                    Real SE = it[slip::n_8c[7]];
                  
                    
                    inn = ( lambda4 * (NW + SE)
                                    + lambda3 * (NE + SW)
                                    + lambda2 * (S + N)
                                    + lambda0 * center)/ g2;
                  }
                else
                  {
                    Real lambda0 = Real(-2.0) * gy2;
                    Real lambda1 = gy2 - gx2;
                    Real lambda3 = Real(0.5) * (gx2 - gxgy);
                    Real lambda4 = Real(0.5) * (gx2 + gxgy);
                   
                    Real center = *it;
                    Real E  = it[slip::n_8c[0]];
                    Real NE = it[slip::n_8c[1]];
                    
                    Real NW = it[slip::n_8c[3]];
                    Real W  = it[slip::n_8c[4]];
                    Real SW = it[slip::n_8c[5]];
                    // Real S  = it[slip::n_8c[6]];
                    Real SE = it[slip::n_8c[7]];
                    inn = ( lambda4 * (NW + SE)
                                  + lambda3 * (NE + SW)
                                  + lambda1 * (E + W)
                                  + lambda0 * center)/ g2;
                  }

      }
 

    return(inn);

  }
};
template<typename Real> struct Inn_Augereau 
{
  Inn_Augereau(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __inn_augereau<Real,_Category>::inn_augereau(it);
  }
};

//
template<typename Real,typename>
struct __inn_alvarez
{
  template <typename _II>
  static Real
  inn_alvarez(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __inn_alvarez<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  inn_alvarez(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
   
    Real inn   = Real(0);

    if(g2 != 0)
      { 

                Real lambda0 = Real(0.25);
                Real lambda1 = Real(2.0) * lambda0 - gx2 / g2;
                Real lambda2 = Real(2.0) * lambda0 - gy2 / g2;
                Real lambda3 = -lambda0 + Real(0.5) * ( Real(1) - gxgy / g2);
                Real lambda4 = -lambda0 + Real(0.5) * ( Real(1) + gxgy / g2);
                lambda0 = Real(-1.0);//-4.0*0.25 
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];
                
                inn = (   lambda4 * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (N + S)
                                + lambda1 * (W + E)
                                + lambda0 * center)/ g2;
      }
  
    return(inn);

  }
};
template<typename Real> struct Inn_Alvarez 
{
  Inn_Alvarez(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __inn_alvarez<Real,_Category>::inn_alvarez(it);
  }
};
//
template<typename Real,typename>
struct __inn_alvarez2
{
  template <typename _II>
  static Real
  inn_alvarez2(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __inn_alvarez2<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  inn_alvarez2(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
   
    Real inn   = Real(0);

    if(g2 != 0)
      { 
                Real absgx = std::abs(gx);
                Real absgy = std::abs(gy);
                Real lambda0 = Real(0);
                if ( (absgx > absgy) && (gxgy>=0) )
                {
                  lambda0 = ( Real(2.0) * gy2 + gx2 + gxgy ) / Real(4.0);
                
                } 
      else if ( (absgx <= absgy) && (gxgy >= 0) )
                {
                  lambda0 = ( Real(2.0) * gx2 + gy2 + gxgy ) / Real(4.0);
                 
                }
      else if ( (absgx > absgy) && (gxgy <= 0) )
                {
                  
                  lambda0 = ( Real(2.0) * gy2 + gx2 + gxgy ) / Real(4.0);       
                 
                }    
      else
                {
                  lambda0 = ( gx2 + Real(2.0) * gy2 + gxgy ) / Real(4.0);
                
                } 

                Real lambda1 = Real(2.0) * lambda0 - gx2;
                Real lambda2  = Real(2.0) *lambda0 - gy2;
                Real lambda3 = -lambda0 + Real(0.5) * ( gx2 - gxgy + gy2 );
                Real lambda4 = -lambda0 + Real(0.5) * ( gx2 + gxgy + gy2 );
                lambda0 = -4.0*lambda0;

                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];
                
                inn = (   lambda4 * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (N + S)
                                + lambda1 * (W + E)
                                + lambda0 * center)/ g2;
      }
  
    return(inn);

  }
};
template<typename Real> struct Inn_Alvarez2 
{
  Inn_Alvarez2(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __inn_alvarez2<Real,_Category>::inn_alvarez2(it);
  }
};
//
template<typename Real,typename>
struct __inn_cohignac
{
  template <typename _II>
  static Real
  inn_cohignac(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __inn_cohignac<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  inn_cohignac(_II it)
  {
    
    Real gx    = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real gy    = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real gx2   = gx * gx;
    Real gy2   = gy * gy;
    Real g2    = gx2 + gy2;
    Real gxgy  = gx * gy;
   
    Real inn   = Real(0);

    if(g2 != 0)
      { 
                Real lambda0 = Real(0.5) - gx2/g2 + (gx2 / g2) * (gx2 / g2);
                Real lambda1 = Real(2.0) * lambda0 - gx2 /g2;
                Real lambda2 = Real(2.0) * lambda0 - gy2 / g2;
                Real lambda3 = -lambda0 + Real(0.5) * (Real(1.0) - gxgy / g2);
                Real lambda4 = -lambda0 + Real(0.5) * (Real(1.0) + gxgy / g2);
                lambda0 = Real(-4.0) * lambda0;
                
                Real center = *it;
                Real E  = it[slip::n_8c[0]];
                Real NE = it[slip::n_8c[1]];
                Real N  = it[slip::n_8c[2]];
                Real NW = it[slip::n_8c[3]];
                Real W  = it[slip::n_8c[4]];
                Real SW = it[slip::n_8c[5]];
                Real S  = it[slip::n_8c[6]];
                Real SE = it[slip::n_8c[7]];
                
                inn = (   lambda4 * (NW + SE)
                                + lambda3 * (NE + SW)
                                + lambda2 * (N + S)
                                + lambda1 * (W + E)
                                + lambda0 * center);
      }
  
    return(inn);

  }
};
template<typename Real> struct Inn_Cohignac
{ 
  Inn_Cohignac(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __inn_cohignac<Real,_Category>::inn_cohignac(it);
  }
};

//
//dvv second derivative in the direction of v
//
template<typename Real,typename>
struct __dvv
{
  template <typename _II, typename Block>
  static Real
  dvv(_II it, Block v)
  {
    assert(Block::SIZE == 1);
    Real dxx = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    if(v[0] < Real(0))
      {
                dxx = -dxx;
      }
    return dxx;
  }
};

template<typename Real>
struct __dvv<Real,std::random_access_iterator_tag>
{
  template <typename _II, typename Block>
   static Real
  dvv(_II it, Block v)
  {
    assert(Block::SIZE == 1);
    return __dvv<Real,std::random_access_iterator_tag>::dvv(it,v);
  }
};

template<typename Real>
struct __dvv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II, typename Block>
   static Real
  dvv(_II it, Block v)
  {
    assert(Block::SIZE == 2);
    Real result = Real(0);
    Real a2 = v[0] * v[0];
    Real b2 = v[1] * v[1];
    Real norm2 = a2 + b2;

    if(norm2 != Real(0))
      {

                Real dxx = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
                Real dyy = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
                Real dxy = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
                
                result = (a2 * dxx + 2.0 * v[0] * v[1] * dxy + b2 * dyy) / norm2 ;
      }
    return result;
  }
};

template<typename Real>
struct __dvv<Real, std::random_access_iterator3d_tag>
{
  template <typename _II, typename Block>
   static Real
  dvv(_II it, Block v)
  {
    assert(Block::SIZE == 3);
    Real result = Real(0);
    Real a2 = v[0] * v[0];
    Real b2 = v[1] * v[1];
    Real c2 = v[2] * v[2];
    Real norm2 = a2 + b2 + c2;
  
    if(norm2 != Real(0))
      {
                Real dxx = __dxx<Real,std::random_access_iterator3d_tag>::dxx(it);
                Real dyy = __dyy<Real,std::random_access_iterator3d_tag>::dyy(it);
                Real dzz = __dzz<Real,std::random_access_iterator3d_tag>::dzz(it);
                
                Real dxy = __dxy<Real,std::random_access_iterator3d_tag>::dxy(it);
                Real dxz = __dxz<Real,std::random_access_iterator3d_tag>::dxz(it);
                Real dyz = __dyz<Real,std::random_access_iterator3d_tag>::dyz(it);

                Real ab  = v[0] * v[1];
                Real ac  = v[0] * v[2];
                Real bc  = v[1] * v[2];
                
                Real a1 = a2 * dxx;
                Real a2 = b2 * dyy;
                Real a3 = c2 * dzz;
                Real a4 = 2.0 * (ab * dxy);
                Real a5 = 2.0 * (ac * dxz);
                Real a6 = 2.0 * (bc * dyz);


                result =  ((a1 + a2) + (a3 + a4) + (a5 + a6)) / norm2;
                
      }
    return result;
  }
};



template<typename Real, typename Block> struct Dvv 
{
  Dvv(Block v)
    :v_(v)
  {}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __dvv<Real,_Category>::dvv(it,v_);
  }
  Block v_;
};


//
template<typename Real,typename>
struct __minmodx
{
  template <typename _II>
  static Real
  minmodx(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
    
    return minmod<Real>(dxp,dxm);
  }
};

template<typename Real>
struct __minmodx<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  minmodx(_II it)
  {
    
    Real dxp = __dx_fwd<Real,std::random_access_iterator2d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator2d_tag>::dx_bck(it);
   return minmod<Real>(dxp,dxm);
    
  }
};

template<typename Real>
struct __minmodx<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  minmodx(_II it)
  {
    
    Real dxp = __dx_fwd<Real,std::random_access_iterator3d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator3d_tag>::dx_bck(it);
   return minmod<Real>(dxp,dxm);
    
  }
};
template<typename Real> struct Minmodx
{ 
  Minmodx(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __minmodx<Real,_Category>::minmodx(it);
  }
};

//
template<typename Real,typename>
struct __minmody
{
  template <typename _II>
  static Real
  minmody(_II it)
  {
    Real dyp = __dy_fwd<Real,std::random_access_iterator_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator_tag>::dy_bck(it);
    
    return minmod<Real>(dyp,dym);
  }
};

template<typename Real>
struct __minmody<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  minmody(_II it)
  {
    
    Real dyp = __dy_fwd<Real,std::random_access_iterator2d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator2d_tag>::dy_bck(it);
   return minmod<Real>(dyp,dym);
    
  }
};

template<typename Real>
struct __minmody<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  minmody(_II it)
  {
    
    Real dyp = __dy_fwd<Real,std::random_access_iterator3d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator3d_tag>::dy_bck(it);
   return minmod<Real>(dyp,dym);
    
  }
};
template<typename Real> struct Minmody
{ 
  Minmody(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __minmody<Real,_Category>::minmody(it);
  }
};


//
template<typename Real,typename>
struct __minmodz
{
  template <typename _II>
  static Real
  minmodz(_II it)
  {
    Real dzp = __dz_fwd<Real,std::random_access_iterator_tag>::dz_fwd(it);
    Real dzm = __dz_bck<Real,std::random_access_iterator_tag>::dz_bck(it);
    
    return minmod<Real>(dzp,dzm);
  }
};

template<typename Real>
struct __minmodz<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  minmodz(_II it)
  {
    
    Real dzp = __dz_fwd<Real,std::random_access_iterator2d_tag>::dz_fwd(it);
    Real dzm = __dz_bck<Real,std::random_access_iterator2d_tag>::dz_bck(it);
   return minmod<Real>(dzp,dzm);
    
  }
};

template<typename Real>
struct __minmodz<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  minmodz(_II it)
  {
    
    Real dzp = __dz_fwd<Real,std::random_access_iterator3d_tag>::dz_fwd(it);
    Real dzm = __dz_bck<Real,std::random_access_iterator3d_tag>::dz_bck(it);
   return minmod<Real>(dzp,dzm);
    
  }
};
template<typename Real> struct Minmodz
{ 
  Minmodz(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __minmodz<Real,_Category>::minmodz(it);
  }
};

//
// grad minus Osher Sethian
//
template<typename Real,typename>
struct __gradminus_OS
{
  template <typename _II>
  static Real
  gradminus_OS(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
    Real max_x = osher_sethian_minus<Real>(dxp,dxm);
    return std::sqrt(max_x);
  }
};

template<typename Real>
struct __gradminus_OS<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  gradminus_OS(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator2d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator2d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator2d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator2d_tag>::dy_bck(it);
    Real max_x = osher_sethian_minus<Real>(dxp,dxm);
    Real max_y = osher_sethian_minus<Real>(dyp,dym);
   return std::sqrt(max_x + max_y);
    
  }
};

template<typename Real>
struct __gradminus_OS<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  gradminus_OS(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator3d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator3d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator3d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator3d_tag>::dy_bck(it);
    Real dzp = __dz_fwd<Real,std::random_access_iterator3d_tag>::dz_fwd(it);
    Real dzm = __dz_bck<Real,std::random_access_iterator3d_tag>::dz_bck(it);
  
    Real max_x = osher_sethian_minus<Real>(dxp,dxm);
    Real max_y = osher_sethian_minus<Real>(dyp,dym);
    Real max_z = osher_sethian_minus<Real>(dzp,dzm);
    
   return std::sqrt(max_x + max_y + max_z);
    
  }
};
template<typename Real> struct Gradminus_OS
{ 
  Gradminus_OS(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __gradminus_OS<Real,_Category>::gradminus_OS(it);
  }
};


//
// grad plus Osher Sethian
//
template<typename Real,typename>
struct __gradplus_OS
{
  template <typename _II>
  static Real
  gradplus_OS(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
    Real max_x = osher_sethian_plus<Real>(dxp,dxm);
    return std::sqrt(max_x);
  }
};

template<typename Real>
struct __gradplus_OS<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  gradplus_OS(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator2d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator2d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator2d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator2d_tag>::dy_bck(it);
    Real max_x = osher_sethian_plus<Real>(dxp,dxm);
    Real max_y = osher_sethian_plus<Real>(dyp,dym);
   return std::sqrt(max_x + max_y);
    
  }
};

template<typename Real>
struct __gradplus_OS<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  gradplus_OS(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator3d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator3d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator3d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator3d_tag>::dy_bck(it);
    Real dzp = __dz_fwd<Real,std::random_access_iterator3d_tag>::dz_fwd(it);
    Real dzm = __dz_bck<Real,std::random_access_iterator3d_tag>::dz_bck(it);
  
    Real max_x = osher_sethian_plus<Real>(dxp,dxm);
    Real max_y = osher_sethian_plus<Real>(dyp,dym);
    Real max_z = osher_sethian_plus<Real>(dzp,dzm);
    
   return std::sqrt(max_x + max_y + max_z);
    
  }
};
template<typename Real> struct Gradplus_OS
{ 
  Gradplus_OS(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __gradplus_OS<Real,_Category>::gradplus_OS(it);
  }
};


//
// grad minus Brockett Maragos
//
template<typename Real,typename>
struct __gradminus_BM
{
  template <typename _II>
  static Real
  gradminus_BM(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
    Real max_x = max_brockett_maragos_minus<Real>(dxp,dxm);
    return std::sqrt(max_x);
  }
};

template<typename Real>
struct __gradminus_BM<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  gradminus_BM(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator2d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator2d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator2d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator2d_tag>::dy_bck(it);
    Real max_x = max_brockett_maragos_minus<Real>(dxp,dxm);
    Real max_y = max_brockett_maragos_minus<Real>(dyp,dym);
   return std::sqrt(max_x + max_y);
    
  }
};

template<typename Real>
struct __gradminus_BM<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  gradminus_BM(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator3d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator3d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator3d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator3d_tag>::dy_bck(it);
    Real dzp = __dz_fwd<Real,std::random_access_iterator3d_tag>::dz_fwd(it);
    Real dzm = __dz_bck<Real,std::random_access_iterator3d_tag>::dz_bck(it);
  
    Real max_x = max_brockett_maragos_minus<Real>(dxp,dxm);
    Real max_y = max_brockett_maragos_minus<Real>(dyp,dym);
    Real max_z = max_brockett_maragos_minus<Real>(dzp,dzm);
    
   return std::sqrt(max_x + max_y + max_z);
    
  }
};
template<typename Real> struct Gradminus_BM
{ 
  Gradminus_BM(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __gradminus_BM<Real,_Category>::gradminus_BM(it);
  }
};

//
// grad plus Brockett Maragos
//
template<typename Real,typename>
struct __gradplus_BM
{
  template <typename _II>
  static Real
  gradplus_BM(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator_tag>::dx_bck(it);
    Real max_x = max_brockett_maragos_plus<Real>(dxp,dxm);
    return std::sqrt(max_x);
  }
};

template<typename Real>
struct __gradplus_BM<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  gradplus_BM(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator2d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator2d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator2d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator2d_tag>::dy_bck(it);
    Real max_x = max_brockett_maragos_plus<Real>(dxp,dxm);
    Real max_y = max_brockett_maragos_plus<Real>(dyp,dym);
   return std::sqrt(max_x + max_y);
    
  }
};

template<typename Real>
struct __gradplus_BM<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  gradplus_BM(_II it)
  {
    Real dxp = __dx_fwd<Real,std::random_access_iterator3d_tag>::dx_fwd(it);
    Real dxm = __dx_bck<Real,std::random_access_iterator3d_tag>::dx_bck(it);
    Real dyp = __dy_fwd<Real,std::random_access_iterator3d_tag>::dy_fwd(it);
    Real dym = __dy_bck<Real,std::random_access_iterator3d_tag>::dy_bck(it);
    Real dzp = __dz_fwd<Real,std::random_access_iterator3d_tag>::dz_fwd(it);
    Real dzm = __dz_bck<Real,std::random_access_iterator3d_tag>::dz_bck(it);
  
    Real max_x = max_brockett_maragos_plus<Real>(dxp,dxm);
    Real max_y = max_brockett_maragos_plus<Real>(dyp,dym);
    Real max_z = max_brockett_maragos_plus<Real>(dzp,dzm);
    
   return std::sqrt(max_x + max_y + max_z);
    
  }
};
template<typename Real> struct Gradplus_BM
{ 
  Gradplus_BM(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __gradplus_BM<Real,_Category>::gradplus_BM(it);
  }
};

//
// grad minus minmod
//
template<typename Real,typename>
struct __gradminus_minmod
{
  template <typename _II>
  static Real
  gradminus_minmod(_II it)
  {
    Real dx = __minmodx<Real,std::random_access_iterator_tag>::minmodx(it);
    return std::sqrt(dx);
  }
};

template<typename Real>
struct __gradminus_minmod<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  gradminus_minmod(_II it)
  {
    Real dx = __minmodx<Real,std::random_access_iterator2d_tag>::minmodx(it);
    Real dy = __minmody<Real,std::random_access_iterator2d_tag>::minmody(it);
    
    return std::sqrt(dx + dy);
  }
};

template<typename Real>
struct __gradminus_minmod<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  gradminus_minmod(_II it)
  {
  
    Real dx = __minmodx<Real,std::random_access_iterator3d_tag>::minmodx(it);
    Real dy = __minmody<Real,std::random_access_iterator3d_tag>::minmody(it);
    Real dz = __minmodz<Real,std::random_access_iterator3d_tag>::minmodz(it);
    
   return std::sqrt(dx + dy + dz);
    
  }
};
template<typename Real> struct Gradminus_minmod
{ 
  Gradminus_minmod(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __gradminus_minmod<Real,_Category>::gradminus_minmod(it);
  }
};

template<typename Real> struct Gradplus_minmod
{ 
  Gradplus_minmod(){}
  template<typename Iterator>
  Real operator() (Iterator it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __gradminus_minmod<Real,_Category>::gradminus_minmod(it);
  }
};


//                   Curvatures
//
// cmin : minimal principal curvature
//
template<typename Real,typename>
struct __min_curv
{
  template <typename _II>
  static Real
  min_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
  }
};

template<typename Real>
struct __min_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  min_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
  }
};

template<typename Real>
struct __min_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  min_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real q = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real p2 = p * p;
    Real q2 = q * q;
    Real H = 1 + p2 + q2;
    Real n = Real(std::sqrt(H));
    Real n3= H * n;
    
    Real n4 = H * H;
    Real ctot = (r * t - s * s) / n4;
    Real rau = ( (Real(1) + p2) * t - Real(2) * p * q * s + (Real(1) + q2) * r) / n3;
    
    return( rau - Real(std::sqrt(rau * rau - ctot)) );
  }
};




template<typename Real> struct Min_Curv 
{
  
  Min_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __min_curv<Real,_Category>::min_curv(it);
  }
};

//
// max_curv : maximal principal curvature
//
template<typename Real,typename>
struct __max_curv
{
  template <typename _II>
  static Real
  max_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
    
  }
};

template<typename Real>
struct __max_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  max_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
  }
};

template<typename Real>
struct __max_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  max_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real q = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real p2 = p * p;
    Real q2 = q * q;
    Real H = 1 + p2 + q2;
    Real n = Real(std::sqrt(H));
    Real n3= H * n;
    
    Real n4 = H * H;
    Real ctot = (r * t - s * s) / n4;
    Real rau = ( (Real(1) + p2) * t - Real(2) * p * q * s + (Real(1) + q2) * r) / n3;
    
    return( rau + Real(std::sqrt(rau * rau - ctot)) );
  }
};





template<typename Real> struct Max_Curv 
{
  
  Max_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __max_curv<Real,_Category>::max_curv(it);
  }
};


//
// mean_curv : mean principal curvature
//
template<typename Real,typename>
struct __mean_curv
{
  template <typename _II>
  static Real
  mean_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
   
  }
};

template<typename Real>
struct __mean_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  mean_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
  }
};

template<typename Real>
struct __mean_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  mean_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real q = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real p2 = p * p;
    Real q2 = q * q;
    Real H = 1 + p2 + q2;
    Real n = Real(std::sqrt(H));
    Real n3= H * n;
        
    return ( ( (Real(1) + p2) * t - Real(2) * p * q  * s + (Real(1) + q2) * r) /(Real(2.0) * n3) );
  }
};

template<typename Real>
struct __mean_curv<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  mean_curv(_II it)
  {
   
    Real Ix  = __dx<Real,std::random_access_iterator3d_tag>::dx(it);
    Real Iy  = __dy<Real,std::random_access_iterator3d_tag>::dy(it);
    Real Iz  = __dz<Real,std::random_access_iterator3d_tag>::dz(it);
    Real Ixy = __dxy<Real,std::random_access_iterator3d_tag>::dxy(it);
    Real Ixz = __dxz<Real,std::random_access_iterator3d_tag>::dxz(it);
    Real Iyz = __dyz<Real,std::random_access_iterator3d_tag>::dyz(it);
    Real Ixx = __dxx<Real,std::random_access_iterator3d_tag>::dxx(it);
    Real Iyy = __dyy<Real,std::random_access_iterator3d_tag>::dyy(it);
    Real Izz = __dzz<Real,std::random_access_iterator3d_tag>::dzz(it);
    
    Real Ix2 = Ix * Ix;
    Real Iy2 = Iy * Iy;
    Real Iz2 = Iz * Iz;
    
    Real IxIy = Ix * Iy;
    Real IxIz = Ix * Iz;
    Real IyIz = Iy * Iz;
 
    Real a = (Real(1.0) + (Ix2 + Iz2));
    Real b = (Real(1.0) + (Iy2 + Iz2));
    Real c = (Real(1.0) + (Ix2 + Iy2));
  
    Real g = (Real(1.0) + Ix2) + (Iy2 + Iz2);
    
    Real den = Real(std::pow(double(g),double(3.0/2.0)));
    
    Real H = ( 
                       (a * Iyy) + (b * Ixx) + (c * Izz) 
                       - Real(2.0) * IxIy * Ixy
                       - Real(2.0) * IxIz * Ixz  
                       - Real(2.0) * IyIz * Iyz ) / den ;

    return (H / Real(3.0));  
  }
};



template<typename Real> struct Mean_Curv 
{
  
  Mean_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __mean_curv<Real,_Category>::mean_curv(it);
  }
};

//
// total_curv : total principal curvature
//
template<typename Real,typename>
struct __total_curv
{
  template <typename _II>
  static Real
  total_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
  }
};

template<typename Real>
struct __total_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  total_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator_tag>::dx(it);
    Real r = __dxx<Real,std::random_access_iterator_tag>::dxx(it);
    Real p2 = p * p;
    return -r / Real(std::sqrt(1.0 + double(p2)));
  }
};

template<typename Real>
struct __total_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  total_curv(_II it)
  {
    Real p  = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real q  = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real s  = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r  = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t  = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real p2 = p * p;
    Real q2 = q * q;

    Real n  = (Real(1) + p2 + q2) * (Real(1) + p2 + q2);
           
    return (( r * t - s * s) / n);
  }
};

template<typename Real>
struct __total_curv<Real, std::random_access_iterator3d_tag>
{
  template <typename _II>
   static Real
  total_curv(_II it)
  {
   
    Real Ix  = __dx<Real,std::random_access_iterator3d_tag>::dx(it);
    Real Iy  = __dy<Real,std::random_access_iterator3d_tag>::dy(it);
    Real Iz  = __dz<Real,std::random_access_iterator3d_tag>::dz(it);
    Real Ixy = __dxy<Real,std::random_access_iterator3d_tag>::dxy(it);
    Real Ixz = __dxz<Real,std::random_access_iterator3d_tag>::dxz(it);
    Real Iyz = __dyz<Real,std::random_access_iterator3d_tag>::dyz(it);
    Real Ixx = __dxx<Real,std::random_access_iterator3d_tag>::dxx(it);
    Real Iyy = __dyy<Real,std::random_access_iterator3d_tag>::dyy(it);
    Real Izz = __dzz<Real,std::random_access_iterator3d_tag>::dzz(it);
    
    Real Ix2 = Ix * Ix;
    Real Iy2 = Iy * Iy;
    Real Iz2 = Iz * Iz;
    
    Real Ixy2 = Ixy * Ixy;
    Real Ixz2 = Ixz * Ixz;
    Real Iyz2 = Iyz * Iyz;

   
    Real g = (Real(1.0) + Ix2) + (Iy2 + Iz2);
  
    Real den = Real(std::pow(double(g),double(2.0/5.0)));
    
    Real K = (  (Ixx * Iyy) * Izz 
                                + (Real(2.0) * Ixy) * (Ixz * Iyz) 
                       - (Ixx * Iyz2) - (Iyy * Ixz2) - (Izz * Ixy2)) / den;

    return K;
  
  }
};



template<typename Real> struct Total_Curv 
{
  
  Total_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __total_curv<Real,_Category>::total_curv(it);
  }
};


//
// lambda1 : first eigen value of the hessian matrix
//
template<typename Real,typename>
struct __lambda1_curv
{
  template <typename _II>
  static Real
  lambda1_curv(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __lambda1_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  lambda1_curv(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __lambda1_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  lambda1_curv(_II it)
  {
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real s2 = s * s;

    return (((r + t) + Real(std::sqrt(double( (r - t) * (r - t)  + 4.0 * s2) ) ) )/ Real(2.0) );
  }
};



template<typename Real> struct Lambda1_Curv 
{
  
  Lambda1_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __lambda1_curv<Real,_Category>::lambda1_curv(it);
  }
};


//
// lambda2 : first eigen value of the hessian matrix
//
template<typename Real,typename>
struct __lambda2_curv
{
  template <typename _II>
  static Real
  lambda2_curv(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __lambda2_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  lambda2_curv(_II it)
  {
    return __dxx<Real,std::random_access_iterator_tag>::dxx(it);
  }
};

template<typename Real>
struct __lambda2_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  lambda2_curv(_II it)
  {
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real s2 = s * s;

    return (((r + t) - Real(std::sqrt(double( (r - t) * (r - t)  + 4.0 * s2) ) ) )/ Real(2.0) );
  }
};




template<typename Real> struct Lambda2_Curv 
{
  
  Lambda2_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __lambda2_curv<Real,_Category>::lambda2_curv(it);
  }
};


//
// k_iso: isophote curvature
//
template<typename Real,typename>
struct __iso_curv
{
  template <typename _II>
  static Real
  iso_curv(_II it)
  {
    return __min_curv<Real,std::random_access_iterator_tag>::min_curv(it);
  }
};

template<typename Real>
struct __iso_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  iso_curv(_II it)
  {
    return __min_curv<Real,std::random_access_iterator_tag>::min_curv(it);
  }
};

template<typename Real>
struct __iso_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  iso_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real q = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real p2 = p * p;
    Real q2 = q * q;
    Real H = p2 + q2;
    Real n = Real(std::sqrt(H));
    Real n3= H * n;
    
    Real result = Real(0.0);
   
    if(n3 != Real(0.0))
      {
                result = (Real(2.0) * p * q * s - q2 * r - p2 * t ) / n3;
      }
    else
      {
                result  = std::numeric_limits<Real>::max();
      }
    return result;
  }
};


template<typename Real> struct Iso_Curv 
{
  
  Iso_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __iso_curv<Real,_Category>::iso_curv(it);
  }
};

//
// courant line curvature
//
template<typename Real,typename>
struct __courant_line_curv
{
  template <typename _II>
  static Real
  courant_line_curv(_II it)
  {
    return __min_curv<Real,std::random_access_iterator_tag>::min_curv(it);
  }
};

template<typename Real>
struct __courant_line_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  courant_line_curv(_II it)
  {
    return __min_curv<Real,std::random_access_iterator_tag>::min_curv(it);
  }
};

template<typename Real>
struct __courant_line_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  courant_line_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real q = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real p2 = p * p;
    Real q2 = q * q;
    Real H = p2 + q2;
    Real n = Real(std::sqrt(H));
    Real n3= H * n;
    
    Real result = Real(0.0);
   
    if(n3 != Real(0.0))
      {
                result = (p * q * (t - r) + s * (p2  - q2) ) / n3;
      }
    else
      {
                result  = std::numeric_limits<Real>::max();
      }
    return result;
  }
};


template<typename Real> struct Courant_Line_Curv 
{
  
  Courant_Line_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __courant_line_curv<Real,_Category>::courant_line_curv(it);
  }
};


//
// computes the relative variation of the gradient along courant lines
//
template<typename Real,typename>
struct __delta_g_curv
{
  template <typename _II>
  static Real
  delta_g_curv(_II it)
  {
    return __min_curv<Real,std::random_access_iterator_tag>::min_curv(it);
  }
};

template<typename Real>
struct __delta_g_curv<Real,std::random_access_iterator_tag>
{
  template <typename _II>
   static Real
  delta_g_curv(_II it)
  {
    return __min_curv<Real,std::random_access_iterator_tag>::min_curv(it);
  }
};

template<typename Real>
struct __delta_g_curv<Real, std::random_access_iterator2d_tag>
{
  template <typename _II>
   static Real
  delta_g_curv(_II it)
  {
    Real p = __dx<Real,std::random_access_iterator2d_tag>::dx(it);
    Real q = __dy<Real,std::random_access_iterator2d_tag>::dy(it);
    Real s = __dxy<Real,std::random_access_iterator2d_tag>::dxy(it);
    Real r = __dxx<Real,std::random_access_iterator2d_tag>::dxx(it);
    Real t = __dyy<Real,std::random_access_iterator2d_tag>::dyy(it);
    Real p2 = p * p;
    Real q2 = q * q;
    Real H = p2 + q2;
    Real n = Real(std::sqrt(H));
    Real n3= H * n;
    
    Real result = Real(0.0);
   
    if(n3 != Real(0.0))
      {
                result = - (((2.0 * p) * (q * s) + (q2 * r) + (p2 * t) ) / n3);
      }
    else
      {
                result  = std::numeric_limits<Real>::max();
      }
    return result;
  }
};


template<typename Real> struct Delta_G_Curv 
{
  
  Delta_G_Curv(){}
  template<typename Iterator>
  Real operator() (Iterator& it) const
  {
    typedef typename std::iterator_traits<Iterator>::iterator_category _Category;
    return __delta_g_curv<Real,_Category>::delta_g_curv(it);
  }
};



/* @{ */
template <typename RandomAccessIteratorNd1, 
                  typename RandomAccessIteratorNd2,
                  typename Operator>
void differential_operator(RandomAccessIteratorNd1 in_first, 
                                                   RandomAccessIteratorNd1 in_last,
                                                   RandomAccessIteratorNd2 out_first,
                                                   Operator& op)
{
  for(; in_first != in_last; ++in_first, ++out_first)
    {
      *out_first = op(in_first);
    }
}
template <typename RandomAccessIteratorNd1, 
                  typename RandomAccessIteratorNd2,
                  typename RandomAccessIteratorNd3,
                  typename Operator1,
                  typename Operator2>
void differential_operator_2(RandomAccessIteratorNd1 in_first, 
                                                     RandomAccessIteratorNd1 in_last,
                                                     RandomAccessIteratorNd2 out_first1,
                                                     RandomAccessIteratorNd3 out_first2,
                                                    
                                                     Operator1& op1,
                                                     Operator2& op2)
{
  for(; in_first != in_last; ++in_first, ++out_first1, ++out_first2)
    {
      *out_first1 = op1(in_first);
      *out_first2 = op2(in_first);
    }
}

template <typename RandomAccessIteratorNd1, 
                  typename RandomAccessIteratorNd2,
                  typename RandomAccessIteratorNd3,
                  typename RandomAccessIteratorNd4,
                  typename Operator1,
                  typename Operator2,
                  typename Operator3>
void differential_operator_3(RandomAccessIteratorNd1 in_first, 
                                                     RandomAccessIteratorNd1 in_last,
                                                     RandomAccessIteratorNd2 out_first1,
                                                     RandomAccessIteratorNd3 out_first2,
                                                     RandomAccessIteratorNd4 out_first3,
                                                     Operator1& op1,
                                                     Operator2& op2,
                                                     Operator3& op3)
{
  for(; in_first != in_last; ++in_first, ++out_first1, ++out_first2, 
++out_first3)
    {
      *out_first1 = op1(in_first);
      *out_first2 = op2(in_first);
      *out_first3 = op3(in_first);
    }
}

/* @} */

}//slip::

#endif //SLIP_DERIVATIVES_HPP