genericpradosschemesforfastmarching_3d.h.svn-base

fast marching method

            }

            // Initialisation of Root.
            T Root[3][3][3];
            for(s1=0; s1<3; s1++) for(s2=0; s2<3; s2++) for(s3=0; s3<3; s3++) {
                Root[s1][s2][s3] = this->big;
            }

            // Initialisation of the optimal dynamics.
            // In dimension 3, the dynamics (function f associated to an HJB eq.) is a 3D-vector.
            // optDymamics1, optDymamics2, optDymamics3 represente the its componantes.
            // As for "Root" and "DoExistSol", we define and compute an optimal
            // Dymamics on each simplex and each "boundaries" of simplexes.

            T optDymamics1[3][3][3];
            T optDymamics2[3][3][3];
            T optDymamics3[3][3][3];

            for(s1=0; s1<3; s1++) for(s2=0; s2<3; s2++) for(s3=0; s3<3; s3++) {
                optDymamics1[s1][s2][s3] = 0;
                optDymamics2[s1][s2][s3] = 0;
                optDymamics3[s1][s2][s3] = 0;
            }


            // ===========================================

            // Computation of Roots such that opt control
            // is in interior of the various Ds1s2s3:

            for(s1=-1; s1<=1; s1+=2)
                for(s2=-1; s2<=1; s2+=2)
                    for(s3=-1; s3<=1; s3+=2) {
                        DoExistSol[s1+1][s2+1][s3+1]
                            = eqSolverOnPart_with_s1s2s3_nonNull(
                                this->_GetValue(x+s1,y,z),
                                this->_GetValue(x,y+s2,z),
                                this->_GetValue(x,y,z+s3),
                                s1,s2,s3,
                                x,y,z,
                                Root[s1+1][s2+1][s3+1],
                                optDymamics1[s1+1][s2+1][s3+1],optDymamics2[s1+1][s2+1][s3+1], optDymamics3[s1+1][s2+1][s3+1]
                                );
                    }

            // ===========================================

            // Computation of Roots such that opt control
            // is in interior of the intersection of Ds1s2s3
            // and set of a such that f_i(x,a) = 0

            s1=0;
            for(s2=-1; s2<=1; s2+=2)
                for(s3=-1; s3<=1; s3+=2) {
                    if ( ! (DoExistSol[-1+1][s2+1][s3+1] || DoExistSol[1+1][s2+1][s3+1]) ) {
                        DoExistSol[s1+1][s2+1][s3+1]
                            = eqSolverOnPart_withOne_si_Null(
                                this->_GetValue(x+s1,y,z),
                                this->_GetValue(x,y+s2,z),
                                this->_GetValue(x,y,z+s3),
                                s1,s2,s3,
                                x,y,z,
                                1,  // Indice_i_Such_siEqualTo0;
                                Root[s1+1][s2+1][s3+1],
                                optDymamics1[s1+1][s2+1][s3+1],optDymamics2[s1+1][s2+1][s3+1], optDymamics3[s1+1][s2+1][s3+1]
                                );
                    }
                    else{
                        DoExistSol[s1+1][s2+1][s3+1] = true;  // useful for the next step!
                    }
                }

            s2=0;
            for(s1=-1; s1<=1; s1+=2)
                for(s3=-1; s3<=1; s3+=2) {
                    if ( ! (DoExistSol[s1+1][-1+1][s3+1] || DoExistSol[s1+1][1+1][s3+1]) ) {
                        DoExistSol[s1+1][s2+1][s3+1]
                            = eqSolverOnPart_withOne_si_Null(
                                this->_GetValue(x+s1,y,z),
                                this->_GetValue(x,y+s2,z),
                                this->_GetValue(x,y,z+s3),
                                s1,s2,s3,
                                x,y,z,
                                2,  // Indice_i_Such_siEqualTo0;
                                Root[s1+1][s2+1][s3+1],
                                optDymamics1[s1+1][s2+1][s3+1],optDymamics2[s1+1][s2+1][s3+1], optDymamics3[s1+1][s2+1][s3+1]
                                );
                    }
                    else{
                        DoExistSol[s1+1][s2+1][s3+1] = true;  // useful for the next step!
                    }
                }

            s3=0;
            for(s2=-1; s2<=1; s2+=2)
                for(s1=-1; s1<=1; s1+=2) {
                    if ( ! (DoExistSol[s1+1][s2+1][-1+1] || DoExistSol[s1+1][s2+1][1+1]) ) {
                        DoExistSol[s1+1][s2+1][s3+1]
                            = eqSolverOnPart_withOne_si_Null(
                                this->_GetValue(x+s1,y,z),
                                this->_GetValue(x,y+s2,z),
                                this->_GetValue(x,y,z+s3),
                                s1,s2,s3,
                                x,y,z,
                                3,  // Indice_i_Such_siEqualTo0;
                                Root[s1+1][s2+1][s3+1],
                                optDymamics1[s1+1][s2+1][s3+1],optDymamics2[s1+1][s2+1][s3+1], optDymamics3[s1+1][s2+1][s3+1]
                                );
                    }
                    else {
                        DoExistSol[s1+1][s2+1][s3+1] = true;  // useful for the next step!
                    }
                }


            // ===========================================
            // computation of Roots such that opt control
            // is in Interior of the intersection of Ds1s2s3
            // and set of a such that f_i(x,a) = 0

            s1=0;
            s2=0;
            for(s3=-1; s3<=1; s3+=2) {
                // More comments could be required here! ....
                if ( ! (
                    DoExistSol[ 0+1][-1+1][s3+1]
                    || DoExistSol[ 0+1][ 1+1][s3+1]
                    || DoExistSol[-1+1][ 0+1][s3+1]
                    || DoExistSol[ 1+1][ 0+1][s3+1]
                    )
                    )
                    {
                        DoExistSol[s1+1][s2+1][s3+1]
                            = eqSolverOnPart_withTwo_si_Null(
                                this->_GetValue(x+s1,y,z),
                                this->_GetValue(x,y+s2,z),
                                this->_GetValue(x,y,z+s3),
                                s1,s2,s3,
                                x,y,z,
                                3,  // Indice_i_Such_siDifferentTo0;
                                Root[s1+1][s2+1][s3+1],
                                optDymamics1[s1+1][s2+1][s3+1],optDymamics2[s1+1][s2+1][s3+1], optDymamics3[s1+1][s2+1][s3+1]
                                );
                    }
            }

            s1=0;
            s3=0;
            for(s2=-1; s2<=1; s2+=2) {
                // More comments could be required here! ....
                if ( ! (
                    DoExistSol[ 0+1][s2+1][-1+1]
                    || DoExistSol[ 0+1][s2+1][ 1+1]
                    || DoExistSol[-1+1][s2+1][ 0+1]
                    || DoExistSol[ 1+1][s2+1][ 0+1]
                    )
                    )
                    {
                        DoExistSol[s1+1][s2+1][s3+1]
                            = eqSolverOnPart_withTwo_si_Null(
                                this->_GetValue(x+s1,y,z),
                                this->_GetValue(x,y+s2,z),
                                this->_GetValue(x,y,z+s3),
                                s1,s2,s3,
                                x,y,z,
                                2,  // Indice_i_Such_siDifferentTo0;
                                Root[s1+1][s2+1][s3+1],
                                optDymamics1[s1+1][s2+1][s3+1],optDymamics2[s1+1][s2+1][s3+1], optDymamics3[s1+1][s2+1][s3+1]
                                );
                    }
            }

            s2=0;
            s3=0;
            for(s1=-1; s1<=1; s1+=2) {
                // More comments could be required here! ....
                if ( ! (
                    DoExistSol[s1+1][ 0+1][-1+1]
                    || DoExistSol[s1+1][ 0+1][ 1+1]
                    || DoExistSol[s1+1][-1+1][ 0+1]
                    || DoExistSol[s1+1][ 1+1][ 0+1]
                    )
                    )
                    {
                        DoExistSol[s1+1][s2+1][s3+1]
                            = eqSolverOnPart_withTwo_si_Null(
                                this->_GetValue(x+s1,y,z),
                                this->_GetValue(x,y+s2,z),
                                this->_GetValue(x,y,z+s3),
                                s1,s2,s3,
                                x,y,z,
                                1,  // Indice_i_Such_siDifferentTo0;
                                Root[s1+1][s2+1][s3+1],
                                optDymamics1[s1+1][s2+1][s3+1],optDymamics2[s1+1][s2+1][s3+1], optDymamics3[s1+1][s2+1][s3+1]
                                );
                    }
            }


            // ===========================================================================


            // The solution is then the smallest Root:
            T minRoot = this->big;      //  the smallest root
            
            for(s1=0; s1<3; s1++)
                for(s2=0; s2<3; s2++)
                    for(s3=0; s3<3; s3++)
                        if ((minRoot > Root[s1][s2][s3]))
                        {
                            minRoot = Root[s1][s2][s3];     // the current min is then Root[s1][s2][s3]
                            gs1=s1; gs2=s2; gs3=s3;         // remember the simplex associated to the current min
                        };


            // Let us record the optimal Dymamics associated to the adequate simplex.
            // The packed optimal dynamics is the one associated to the used simplex.
            // Note: Preservation 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 optimal Dymamics is usefull in
            // many application as for example fibers tracking in DTI (see papers
            // associated to this code).

            setOptDynamics(x,y,z, optDymamics1[gs1][gs2][gs3],optDymamics2[gs1][gs2][gs3],optDymamics3[gs1][gs2][gs3]);

            // we can then return the value for update...
            return minRoot;

        }// End of the methode "_UpdateValue".
    };  // End of the definition of the class.

} // End of namespace FastLevelSet.

#endif