/*------------------------------------------------------------------------------------------------------    File        : GenericPradosSchemesForFastMarching_3D.h   (GCM Library)  Authors      : Emmanuel Prados (UCLA/INRIA), Christophe Lenglet (INRIA), Jean-Philippe Pons (INRIA)  Description : This method explicits the scheme proposed by Prados etal. in dimension 3 (this scheme is                shortly described in sections 4.3 of INRIA Research Report 5845  -- in particular, see 		sections 4.3.1 and 4.3.2).    --------------  License     : This software is governed by the CeCILL-C license under French law and abiding by the   rules of distribution of free software.   Users can use, modify and/ or redistribute the software under the terms of the CeCILL-C. In particular,   the exercising of this right is conditional upon the obligation to make available to the community the   modifications made to the source code of the software so as to contribute to its evolution (e.g. by the   mean of the web; i.e. by publishing a web page).    In this respect, the risks associated with loading, using, modifying and/or developing or reproducing   the software by the user are brought to the user's attention, given its Free Software status, which may   make it complicated to use, with the result that its use is reserved for developers and experienced   professionals having in-depth computer knowledge. Users are therefore encouraged to load and test the   suitability of the software 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 of   security. This Agreement may be freely reproduced and published, provided it is not altered, and that   no provisions are either added or removed herefrom.     CeCILL-C FREE SOFTWARE LICENSE AGREEMENT is available in the file                           Licence_CeCILL-C_V1-en.txt   or at                            http://www.cecill.info/index.en.html.  This Agreement may apply to any or all software for which the holder of the economic rights decides to   submit the use thereof to its provisions.    --------------  Associated publications  : This C++ code corresponds to the implementation of the algorithm presented   in the following articles:    - E. Prados, C. Lenglet, J.P. Pons, N. Wotawa, R. Deriche, O. Faugeras, S. Soatto;     Control Theory and Fast Marching Methods for Brain Connectivity Mapping; INRIA Research Report 5845     UCLA Computer Science Department Technical Report 060004, February 2006.  - E. Prados, C. Lenglet, J.P. Pons, N. Wotawa, R. Deriche, O. Faugeras, S. Soatto;     Control Theory and Fast Marching Methods for Brain Connectivity Mapping;     Proc. IEEE Computer Society Conference on Computer Vision and Pattern Recognition, New York, NY,     I: 1076-1083, June 17-22, 2006.    - For more references, we refer to the official web page of the GCM Library and to authors' web pages.  Please, if you use the GCM library in you work, make sure you will include the reference to the work   of the authors in your publications.    ----------------------------------------------------------------------------------------------------*/#ifndef GENERICPRADOSSCHEMESFORFASTMARCHING_3D_H#define GENERICPRADOSSCHEMESFORFASTMARCHING_3D_H#include "Globals.h"#include "FastMarching_WithOptimalDynamics.h"/**************************************************/// Dimension 3 =====================================/**************************************************//*Here, we deal with the equation of the form:$$ H(x,\nabla u(x)) = 0 . $$Let $H^*$ be the Legendre transform of $H$.For $i \in [1..N]$, we denote$$ H_i(x,p) = sup_{a\in Dom(H^*), a_i = 0}  a \cdot p - H^*(x,a)$$.$$ H_ij(x,p) = sup_{a\in Dom(H^*), a_i = 0 and a_j = 0,  j \neq i}  a \cdot p - H^*(x,a)$$.*//**************************************************/// Note about the optimal dynamics: // The preservation and the transmission of the optimal dymamics is not usefull if we only want // to compute the viscosity solution of the considered equation by the Fast Marching Method. // Nevertheless the computation of the optimal dymamics (function $f$ in Prados's papers) is // necessary inside (i.e. when we develop) the functions//      bool eqSolverOnPart_with_s1s2s3_nonNull(...)//      bool eqSolverOnPart_withOne_si_Null(...)//      bool eqSolverOnPart_withTwo_si_Null(...)// in order to be able to know if the considered simplex is a good candidate.// Note: preservation and transmission of the optimal dymamics can be usefull in many applications// as for example fibers tracking in DTI.namespace FastLevelSet {    template <typename T = float>    class PradosSchemesForFastMarching_3D : public FastMarching_WithOptDynamics<T> {    public:        // Constructor        PradosSchemesForFastMarching_3D(T *_data, int _width, int _height, int _depth, double *_voro = NULL) : FastMarching_WithOptDynamics<T>(_data,_width,_height,_depth,_voro) {}        // Destructor        virtual ~PradosSchemesForFastMarching_3D() {}    protected:        //////////////////////////////////////////////////////        // These next THREE Functions must be overcharged !        //////////////////////////////////////////////////////        //////////////////////////////////////////////////////////////////////////        //        //  Let us denote:        //        //      $$ H_{s1s2s3}(x,p) = sup_{a\in Dom_{s1s2s3}}  a \cdot p - H^*(x,a)$$.        //      where        //      $$ Dom_{s1s2s3}  = \{ a\in Dom(H^*) such that for all i, sign(a_i)*si < 0,        //        //  eqSolverOnPart_with_s1s2s3_nonNull must solve the following equation (in $t$) :        //        //  $ H_{s1s2s3}(x,p_t) = 0 $     (equation (1)),        //        //  where $[p_t]_i = [ t-u(x+sihiei)] / (-sihi)$.        //        //  Solve equation (1) is equivalent to solve equation        //        //  $ H(x,p_t) = 0 $     (equation (2)),        //        //  and then, amongst the solutions of equation (2),  chose        //  the solution $t_0$ of such that:        //  $\nabla H (x,p_{t_0}) * si <0 $.        //        //  Notes:        //  1)  $\nabla H (x,p_{t_0})$ corresponds with the optimal        //      control of $ H(x,p_{t_0}) = 0 $.        //        //  2)  Here, the $hi$ are not given as parameters, since they        //  are supposed constant, and known in the next inherited        //  classes...        //        //  3) the parameter Ui corresponds with the value $u(x+sihiei)$        //        //////////////////////////////////////////////////////////////////////////        virtual bool eqSolverOnPart_with_s1s2s3_nonNull(            const T U1, const T U2, const T U3,                 // Values of the solution at the considered neigborhood voxels,            const int s1, const int s2, const int s3,           // signs associated to the considered sector,            const int x, const int y, const int z,              // coordinates of the considered voxel,            T &Root,                                            // solution.            T &optDymamics1, T &optDymamics2, T &optDymamics3   // optimal dynamic associated to the solution.            ) const = 0;        virtual bool eqSolverOnPart_withOne_si_Null(            const T U1, const T U2, const T U3,            const int s1, const int s2, const int s3,            const int x, const int y, const int z,            const int indice_si_EqualZero,            T &Root,            T &optDymamics1,    T &optDymamics2,    T &optDymamics3            ) const = 0;        virtual bool eqSolverOnPart_withTwo_si_Null(            const T U1, const T U2, const T U3,            const int s1, const int s2, const int s3,            const int x, const int y, const int z,            const int indice_si_DiffZero,            T &Root,            T &optDymamics1,    T &optDymamics2,    T &optDymamics3            ) const = 0;        //////////////////////////////////////////////////////        // The next function must NOT be overcharged !        //////////////////////////////////////////////////////        ////////////////////////////////////////////////////////////////////////////////////////////////////////////        // Update of a point        ////////////////////////////////////////////////////////////////////////////////////////////////////////////        virtual T _UpdateValue(const int x, const int y, const int z) const        {            // Let us deal with the boundary:            if ((x==0) || (x==this->width-1)  ||                (y==0) || (y==this->height-1) ||                (z==0) || (z==this->depth-1))                return this->big;            int gs1=0, gs2=0, gs3=0;    //  give the simplex of the smallest                                        //  root (usefull for recovering the optimal dynamics).            return _UpdateValueB(x,y,z,gs1,gs2,gs3);                    }// End of the methode "_UpdateValue".                // With Voronoi information        virtual T _UpdateValue(const int x, const int y, const int z, int &mx, int& my, int &mz) const        {            // Let us deal with the boundary:            if ((x==0) || (x==this->width-1)  ||                (y==0) || (y==this->height-1) ||                (z==0) || (z==this->depth-1))                return this->big;                        int gs1=0, gs2=0, gs3=0;    //  give the simplex of the smallest                                        //  root (usefull for recovering the optimal dynamics).            T val = _UpdateValueB(x,y,z,gs1,gs2,gs3);            // Voronoi information            mx = x + gs1-1;            my = y;            mz = z;            T res = this->_GetValue(mx,y,z);            if (this->_GetValue(x,gs2+y-1,z) < res)            {                mx = x;                my = y + gs2-1;            }            if (this->_GetValue(x,y,gs3+z-1) < res)            {                my = y;                mz = z + gs3-1;            }            return val;        }// End of the methode "_UpdateValue".                virtual T _UpdateValueB(const int x, const int y, const int z, int &gs1, int& gs2, int &gs3) const        {            // Initialisation of s1,s2,S3.            int s1=0, s2=0, s3=0;            // Initialisation of DoExistSol.            bool DoExistSol[3][3][3];            for(s1=0; s1<3; s1++) for(s2=0; s2<3; s2++) for(s3=0; s3<3; s3++) {                DoExistSol[s1][s2][s3]=false;