"Fossies" - the Fresh Open Source Software Archive

Member "mathmod-branches-r508-trunk/pariso/isosurface/Iso3D.cpp" (8 Mar 2021, 132501 Bytes) of package /linux/misc/mathmod-11.0-source.zip:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Iso3D.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 10.1_vs_11.0.

    1 /***************************************************************************
    2 *   Copyright (C) 2021 by Abderrahman Taha                                *
    3 *                                                                         *
    4 *                                                                         *
    5 *   This program is free software; you can redistribute it and/or modify  *
    6 *   it under the terms of the GNU General Public License as published by  *
    7 *   the Free Software Foundation; either version 2 of the License, or     *
    8 *   (at your option) any later version.                                   *
    9 *                                                                         *
   10 *   This program is distributed in the hope that it will be useful,       *
   11 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
   12 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
   13 *   GNU General Public License for more details.                          *
   14 *                                                                         *
   15 *   You should have received a copy of the GNU General Public License     *
   16 *   along with this program; if not, write to the                         *
   17 *   Free Software Foundation, Inc.,                                       *
   18 *   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA            *
   19 ***************************************************************************/
   20 #include "TableMap.h"
   21 #include "Iso3D.h"
   22 #include "internalfunctions.cpp"
   23 #include <QElapsedTimer>
   24 
   25 static uint NbPolyMin;
   26 static Voxel *GridVoxelVarPt;
   27 static double *Results;
   28 static uint NbVertexTmp = 0;
   29 static std::vector<float> NormOriginaltmpVector;
   30 uint NbTriangleIsoSurface,NbPointIsoMap;
   31 uint OrignbX, OrignbY, OrignbZ;
   32 uint Stack_Factor=OrignbX*OrignbY*OrignbZ;
   33 std::vector<float> NormVertexTabVector;
   34 std::vector<uint>  IndexPolyTabMinVector;
   35 std::vector<uint>  IndexPolyTabMinVector2;
   36 std::vector<uint>  IndexPolyTabVector;
   37 static CellNoise *NoiseFunction = new CellNoise();
   38 static ImprovedNoise *PNoise = new ImprovedNoise(4., 4., 4.);
   39 static QElapsedTimer times;
   40 static double IsoComponentId=0;
   41 static int nbvariables=0;
   42 
   43 double CurrentIsoCmpId(const double* p)
   44 {
   45     return((int (p[0]))== 0 ? IsoComponentId:0);
   46 }
   47 
   48 extern double TurbulenceWorley(const double *p)
   49 {
   50     return double (
   51                NoiseFunction->CellNoiseFunc(
   52                    float (p[0]),
   53                    float (p[1]),
   54                    float (p[2]),
   55                    int (p[3]),
   56                    int (p[4]),
   57                    int (p[5]))
   58            );
   59 }
   60 double TurbulencePerlin(const double *p)
   61 {
   62     return double (
   63                PNoise->FractalNoise3D(
   64                    float (p[0]),
   65                    float (p[1]),
   66                    float (p[2]),
   67                    int (p[3]),
   68                    float (p[4]),
   69                    float (p[5])));
   70 }
   71 double MarblePerlin(const double *p)
   72 {
   73     return double (
   74                PNoise->Marble(
   75                    float (p[0]),
   76                    float (p[1]),
   77                    float (p[2]),
   78                    int (p[3])));
   79 }
   80 void IsoWorkerThread::run()
   81 {
   82     IsoCompute(CurrentComponent);
   83 }
   84 void Iso3D::run()
   85 {
   86     BuildIso();
   87 }
   88 void IsoWorkerThread::emitMySignal()
   89 {
   90     emit mySignal(signalVal);
   91 }
   92 void Iso3D::emitErrorSignal()
   93 {
   94     emit ErrorSignal(int(messageerror));
   95 }
   96 void Iso3D::emitUpdateMessageSignal()
   97 {
   98     emit UpdateMessageSignal(message);
   99 }
  100 void Iso3D::BuildIso()
  101 {
  102     IsoBuild(
  103         &(localScene->ArrayNorVer_localPt),
  104         &(localScene->PolyIndices_localPt),
  105         &localScene->PolyNumber,
  106         &(localScene->VertxNumber),
  107         &(localScene->PolyIndices_localPtMin),
  108         &(localScene->NbPolygnNbVertexPtMin),
  109         &(localScene->NbPolygnNbVertexPtMinSize),
  110         &(localScene->componentsinfos));
  111 }
  112 IsoMasterThread::IsoMasterThread()
  113 {
  114     ParisoCondition = -1;
  115     ImplicitFunction =  "cos(x) + cos(y) + cos(z)";
  116     XlimitSup = "4";
  117     YlimitSup = "4";
  118     ZlimitSup = "4";
  119     XlimitInf   = "-4";
  120     YlimitInf   = "-4";
  121     ZlimitInf   = "-4";
  122     Lacunarity = 0.5;
  123     Gain = 1.0;
  124     Octaves = 4;
  125     Cstparser.AddConstant("pi", PI);
  126     componentsNumber = ConstSize = FunctSize = RgbtSize = VRgbtSize = 0;
  127 }
  128 void IsoMasterThread::IsoMasterTable()
  129 {
  130     UsedFunct    = new bool[0];
  131     UsedFunct2   = new bool[0];
  132 }
  133 IsoMasterThread::~IsoMasterThread()
  134 {
  135     if(ParsersAllocated)
  136     {
  137         delete[] implicitFunctionParser;
  138         delete[] xSupParser;
  139         delete[] ySupParser;
  140         delete[] zSupParser;
  141         delete[] xInfParser;
  142         delete[] yInfParser;
  143         delete[] zInfParser;
  144         delete[] Fct;
  145         delete[] ParisoConditionParser;
  146         delete[] VRgbtParser;
  147         delete[] RgbtParser;
  148         delete GradientParser;
  149         delete NoiseParser;
  150         ParsersAllocated = false;
  151     }
  152     delete[] UsedFunct;
  153     delete[] UsedFunct2;
  154     grid.clear();
  155     SliderValues.clear();
  156     SliderNames.clear();
  157     Consts.clear();
  158     ConstNames.clear();
  159     ConstValues.clear();
  160     Rgbts.clear();
  161     RgbtNames.clear();
  162     VRgbts.clear();
  163     VRgbtNames.clear();
  164     Functs.clear();
  165     FunctNames.clear();
  166 }
  167 IsoWorkerThread::~IsoWorkerThread()
  168 {
  169     if(ParsersAllocated)
  170     {
  171         delete[] implicitFunctionParser;
  172         delete[] Fct;
  173         ParsersAllocated = false;
  174     }
  175 }
  176 IsoWorkerThread::IsoWorkerThread()
  177 {
  178     XYZgrid = 40;
  179     stepMorph = 0;
  180     pace = 1.0/30.0;
  181     activeMorph = -1;
  182     AllComponentTraited = false;
  183     ParsersAllocated = false;
  184     StopCalculations = false;
  185 }
  186 void Iso3D::WorkerThreadCopy(IsoWorkerThread *WorkerThreadsTmp)
  187 {
  188     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  189     {
  190         WorkerThreadsTmp[nbthreads].XYZgrid = masterthread->XYZgrid;
  191         WorkerThreadsTmp[nbthreads].MyIndex  = nbthreads+1;
  192         WorkerThreadsTmp[nbthreads].WorkerThreadsNumber = WorkerThreadsNumber;
  193         WorkerThreadsTmp[nbthreads].GridVal = masterthread->GridVal;
  194     }
  195 }
  196 void Iso3D::UpdateThredsNumber(uint NewThreadsNumber)
  197 {
  198     uint OldWorkerThreadsNumber = WorkerThreadsNumber;
  199     WorkerThreadsNumber = ThreadsNumber = NewThreadsNumber;
  200     IsoWorkerThread *workerthreadstmp = new IsoWorkerThread[WorkerThreadsNumber-1];
  201     WorkerThreadCopy(workerthreadstmp);
  202     //Free old memory:
  203     for(uint i=0; i+1<OldWorkerThreadsNumber; i++)
  204         workerthreads[i].DeleteWorkerParsers();
  205     if(OldWorkerThreadsNumber >1)
  206         delete[] workerthreads;
  207     //Assigne newly allocated memory
  208     workerthreads = workerthreadstmp;
  209     masterthread->WorkerThreadsNumber = WorkerThreadsNumber;
  210 }
  211 ErrorMessage  Iso3D::IsoMorph()
  212 {
  213     ErrorMessage err = masterthread->ParserIso();
  214     if(err.iErrorIndex < 0)
  215         ThreadParsersCopy();
  216     return err;
  217 }
  218 ErrorMessage Iso3D::ThreadParsersCopy()
  219 {
  220     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  221     {
  222         workerthreads[nbthreads].xLocal2.clear();
  223         workerthreads[nbthreads].xLocal2 = masterthread->xLocal2;
  224         workerthreads[nbthreads].yLocal2.clear();
  225         workerthreads[nbthreads].yLocal2 = masterthread->yLocal2;
  226         workerthreads[nbthreads].zLocal2.clear();
  227         workerthreads[nbthreads].zLocal2 = masterthread->zLocal2;
  228         workerthreads[nbthreads].activeMorph = masterthread->activeMorph;
  229         workerthreads[nbthreads].AllComponentTraited = masterthread->AllComponentTraited;
  230         workerthreads[nbthreads].XYZgrid = masterthread->XYZgrid;
  231         workerthreads[nbthreads].GridVal = masterthread->GridVal;
  232     }
  233     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  234         workerthreads[nbthreads].DeleteWorkerParsers();
  235     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  236         workerthreads[nbthreads].AllocateParsersForWorkerThread(masterthread->componentsNumber, masterthread->FunctSize);
  237     return(parse_expression2());
  238 }
  239 ErrorMessage  Iso3D::parse_expression2()
  240 {
  241     ErrorMessage NodError;
  242     // Functions :
  243     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  244     {
  245         for(uint ij=0; ij<masterthread->FunctSize; ij++)
  246         {
  247             workerthreads[nbthreads].Fct[ij].AddFunction("CmpId",CurrentIsoCmpId, 1);
  248             workerthreads[nbthreads].Fct[ij].AddConstant("pi", PI);
  249             workerthreads[nbthreads].Fct[ij].AddConstant("ThreadId", workerthreads[nbthreads].MyIndex);
  250         }
  251         for(uint ii=0; ii<masterthread->FunctSize; ii++)
  252         {
  253             for(uint jj=0; jj<masterthread->ConstSize; jj++)
  254             {
  255                 workerthreads[nbthreads].Fct[ii].AddConstant(masterthread->ConstNames[jj], masterthread->ConstValues[jj]);
  256             }
  257             //Add predefined constatnts:
  258             for(uint kk=0; kk<masterthread->Nb_Sliders; kk++)
  259             {
  260                 workerthreads[nbthreads].Fct[ii].AddConstant(masterthread->SliderNames[kk], masterthread->SliderValues[kk]);
  261             }
  262         }
  263         for(uint ii=0; ii<masterthread->FunctSize; ii++)
  264         {
  265             workerthreads[nbthreads].Fct[ii].AddFunction("NoiseW",TurbulenceWorley, 6);
  266             workerthreads[nbthreads].Fct[ii].AddFunction("fhelix1",fhelix1, 10);
  267             workerthreads[nbthreads].Fct[ii].AddFunction("fhelix2",fhelix2, 10);
  268             workerthreads[nbthreads].Fct[ii].AddFunction("f_hex_y",f_hex_y, 4);
  269             workerthreads[nbthreads].Fct[ii].AddFunction("p_skeletal_int",p_skeletal_int, 3);
  270             workerthreads[nbthreads].Fct[ii].AddFunction("fmesh",fmesh, 8);
  271             workerthreads[nbthreads].Fct[ii].AddFunction("NoiseP",TurbulencePerlin, 6);
  272             workerthreads[nbthreads].Fct[ii].AddFunction("MarbleP",MarblePerlin, 4);
  273             for(uint jj=0; jj<ii; jj++)
  274                 if(masterthread->UsedFunct2[ii*masterthread->FunctSize+jj])
  275                     workerthreads[nbthreads].Fct[ii].AddFunction(masterthread->FunctNames[jj], workerthreads[nbthreads].Fct[jj]);
  276             if ((masterthread->stdError.iErrorIndex = workerthreads[nbthreads].Fct[ii].Parse(masterthread->Functs[ii],"x,y,z,t")) >= 0)
  277             {
  278                 masterthread->stdError.strError = masterthread->Functs[ii];
  279                 return masterthread->stdError;
  280             }
  281             workerthreads[nbthreads].Fct[ii].AllocateStackMemory(Stack_Factor, nbvariables);
  282         }
  283     }
  284     //Add defined constantes:
  285     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  286     {
  287         for(uint i=0; i<masterthread->componentsNumber; i++)
  288         {
  289             workerthreads[nbthreads].implicitFunctionParser[i].AddConstant("pi", PI);
  290             workerthreads[nbthreads].implicitFunctionParser[i].AddConstant("ThreadId", workerthreads[nbthreads].MyIndex);
  291             for(uint j=0; j<masterthread->ConstSize; j++)
  292             {
  293                 workerthreads[nbthreads].implicitFunctionParser[i].AddConstant(masterthread->ConstNames[j], masterthread->ConstValues[j]);
  294             }
  295             //Add predefined sliders constatnts:
  296             for(uint k=0; k<masterthread->Nb_Sliders; k++)
  297             {
  298                 workerthreads[nbthreads].implicitFunctionParser[i].AddConstant(masterthread->SliderNames[k], masterthread->SliderValues[k]);
  299             }
  300         }
  301     }
  302     // Add defined functions :
  303     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  304     {
  305         for(uint i=0; i<masterthread->componentsNumber; i++)
  306         {
  307             //Functions:
  308             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
  309             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("fhelix1",fhelix1, 10);
  310             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("fhelix2",fhelix2, 10);
  311             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("f_hex_y",f_hex_y, 4);
  312             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("p_skeletal_int",p_skeletal_int, 3);
  313             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("fmesh",fmesh, 8);
  314             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
  315             workerthreads[nbthreads].implicitFunctionParser[i].AddFunction("MarbleP",MarblePerlin, 4);
  316 
  317             for(uint j=0; j<masterthread->FunctSize; j++)
  318             {
  319                 if(masterthread->UsedFunct[i*masterthread->FunctSize+j])
  320                 {
  321                     workerthreads[nbthreads].implicitFunctionParser[i].AddFunction(masterthread->FunctNames[j], workerthreads[nbthreads].Fct[j]);
  322                 }
  323             }
  324         }
  325     }
  326     // Final step: parsing
  327     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  328     {
  329         for(uint index=0; index< masterthread->componentsNumber; index++)
  330         {
  331             if ((masterthread->stdError.iErrorIndex = workerthreads[nbthreads].implicitFunctionParser[index].Parse(masterthread-> ImplicitStructs[index].fxyz, "x,y,z,t,i_indx,j_indx,k_indx,max_ijk")) >= 0)
  332             {
  333                 masterthread->stdError.strError = masterthread->ImplicitStructs[index].fxyz;
  334                 return masterthread->stdError;
  335             }
  336         }
  337     }
  338     return NodError;
  339 }
  340 Iso3D::Iso3D( uint nbThreads,
  341               uint nbGrid,
  342               uint factX,
  343               uint factY,
  344               uint factZ)
  345 {
  346     OrignbX= factX;
  347     OrignbY= factY;
  348     OrignbZ=factZ;
  349     Stack_Factor = factX*factY*factZ;
  350     NbTriangleIsoSurface = 0;
  351     NbPointIsoMap = 0;
  352     WorkerThreadsNumber = ThreadsNumber = nbThreads;
  353     masterthread  = new IsoMasterThread();
  354     masterthread->IsoMasterTable();
  355     workerthreads = new IsoWorkerThread[WorkerThreadsNumber-1];
  356     masterthread->XYZgrid = nbGrid;
  357     masterthread->GridVal = nbGrid;
  358     masterthread->MyIndex = 0;
  359     masterthread->WorkerThreadsNumber = WorkerThreadsNumber;
  360     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  361     {
  362         workerthreads[nbthreads].XYZgrid = nbGrid;
  363         workerthreads[nbthreads].GridVal = nbGrid;
  364         workerthreads[nbthreads].MyIndex = nbthreads+1;
  365         workerthreads[nbthreads].WorkerThreadsNumber = WorkerThreadsNumber;
  366     }
  367 }
  368 uint IsoMasterThread::HowManyVariables(std::string NewVariables, uint type)
  369 {
  370     std::string tmp, tmp2,tmp3;
  371     size_t position =0, jpos;
  372     uint Nb_variables =0;
  373     while( NewVariables!= "")
  374     {
  375         if((position = NewVariables.find(";")) != std::string::npos)
  376         {
  377             tmp = NewVariables;
  378             tmp2= tmp3 = (tmp.substr(0,position));
  379             jpos = tmp2.find("=");
  380             if(type == 1)
  381             {
  382                 ConstNames.push_back(tmp2.substr(0,jpos));
  383                 Consts.push_back(tmp3.substr(jpos+1,position-1));
  384             }
  385             else if(type == 2)
  386             {
  387                 FunctNames.push_back(tmp2.substr(0,jpos));
  388                 Functs.push_back(tmp3.substr(jpos+1,position-1));
  389             }
  390             else if(type == 3)
  391             {
  392                 RgbtNames.push_back(tmp2.substr(0,jpos));
  393                 Rgbts.push_back(tmp3.substr(jpos+1,position-1));
  394             }
  395             else if(type == 4)
  396             {
  397                 VRgbtNames.push_back(tmp2.substr(0,jpos));
  398                 VRgbts.push_back(tmp3.substr(jpos+1,position-1));
  399             }
  400             tmp2 = NewVariables.substr(position+1, NewVariables.length()-1);
  401             NewVariables = tmp2;
  402             Nb_variables++;
  403         }
  404         else
  405         {
  406             tmp = tmp2 = tmp3 = NewVariables;
  407             jpos = tmp2.find("=");
  408             if(type == 1)
  409             {
  410                 ConstNames.push_back(tmp2.substr(0, jpos));
  411                 Consts.push_back(tmp3.substr(jpos+1,position-1));
  412             }
  413             else if(type == 2)
  414             {
  415                 FunctNames.push_back(tmp2.substr(0, jpos));
  416                 Functs.push_back(tmp3.substr(jpos+1,position-1));
  417             }
  418             else if(type == 3)
  419             {
  420                 RgbtNames.push_back(tmp2.substr(0, jpos));
  421                 Rgbts.push_back(tmp3.substr(jpos+1,position-1));
  422             }
  423             else if(type == 4)
  424             {
  425                 VRgbtNames.push_back(tmp2.substr(0, jpos));
  426                 VRgbts.push_back(tmp3.substr(jpos+1,position-1));
  427             }
  428             NewVariables = "";
  429             Nb_variables++;
  430         }
  431     }
  432     return Nb_variables;
  433 }
  434 uint IsoMasterThread::HowManyIsosurface(std::string ImplicitFct, uint type)
  435 {
  436     std::string tmp, tmp2;
  437     size_t position =0;
  438     uint Nb_implicitfunction =0;
  439     if(type ==0)
  440     {
  441         ImplicitStructs[0].fxyz = ImplicitFct;
  442         while( ImplicitFct!= "")
  443         {
  444             if((position = ImplicitFct.find(";")) !=std::string::npos   )
  445             {
  446                 tmp = ImplicitFct;
  447                 ImplicitStructs[Nb_implicitfunction].fxyz = (tmp.substr(0,position));
  448                 Nb_implicitfunction++;
  449                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  450                 ImplicitFct = tmp2;
  451             }
  452             else
  453             {
  454                 ImplicitStructs[Nb_implicitfunction].fxyz  = ImplicitFct;
  455                 ImplicitFct ="";
  456             }
  457         }
  458         return Nb_implicitfunction;
  459     }
  460     else if(type == 1)  //xmin
  461     {
  462         ImplicitStructs[0].xmin = ImplicitFct;
  463         while( ImplicitFct!= "")
  464         {
  465             if((position = ImplicitFct.find(";")) != std::string::npos)
  466             {
  467                 tmp = ImplicitFct;
  468                 ImplicitStructs[Nb_implicitfunction].xmin = (tmp.substr(0,position));
  469                 Nb_implicitfunction++;
  470                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  471                 ImplicitFct = tmp2;
  472             }
  473             else
  474             {
  475                 ImplicitStructs[Nb_implicitfunction].xmin  = ImplicitFct;
  476                 ImplicitFct ="";
  477             }
  478         }
  479         return Nb_implicitfunction;
  480     }
  481     else if(type == 2)  //xmax
  482     {
  483         ImplicitStructs[0].xmax = ImplicitFct;
  484         while( ImplicitFct!= "")
  485         {
  486             if((position = ImplicitFct.find(";")) != std::string::npos)
  487             {
  488                 tmp = ImplicitFct;
  489                 ImplicitStructs[Nb_implicitfunction].xmax = (tmp.substr(0,position));
  490                 Nb_implicitfunction++;
  491                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  492                 ImplicitFct = tmp2;
  493             }
  494             else
  495             {
  496                 ImplicitStructs[Nb_implicitfunction].xmax  = ImplicitFct;
  497                 ImplicitFct ="";
  498             }
  499         }
  500         return Nb_implicitfunction;
  501     }
  502     else if(type == 3) //ymin
  503     {
  504         ImplicitStructs[0].ymin = ImplicitFct;
  505         while( ImplicitFct!= "")
  506         {
  507             if((position = ImplicitFct.find(";")) != std::string::npos)
  508             {
  509                 tmp = ImplicitFct;
  510                 ImplicitStructs[Nb_implicitfunction].ymin = (tmp.substr(0,position));
  511                 Nb_implicitfunction++;
  512                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  513                 ImplicitFct = tmp2;
  514             }
  515             else
  516             {
  517                 ImplicitStructs[Nb_implicitfunction].ymin  = ImplicitFct;
  518                 ImplicitFct ="";
  519             }
  520         }
  521         return Nb_implicitfunction;
  522     }
  523     else if(type == 4) //ymax
  524     {
  525         ImplicitStructs[0].ymax = ImplicitFct;
  526         while( ImplicitFct!= "")
  527         {
  528             if((position = ImplicitFct.find(";")) != std::string::npos)
  529             {
  530                 tmp = ImplicitFct;
  531                 ImplicitStructs[Nb_implicitfunction].ymax = (tmp.substr(0,position));
  532                 Nb_implicitfunction++;
  533                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  534                 ImplicitFct = tmp2;
  535             }
  536             else
  537             {
  538                 ImplicitStructs[Nb_implicitfunction].ymax  = ImplicitFct;
  539                 ImplicitFct ="";
  540             }
  541         }
  542         return Nb_implicitfunction;
  543     }
  544     else if(type == 5) //zmin
  545     {
  546         ImplicitStructs[0].zmin = ImplicitFct;
  547         while( ImplicitFct!= "")
  548         {
  549             if((position = ImplicitFct.find(";")) != std::string::npos)
  550             {
  551                 tmp = ImplicitFct;
  552                 ImplicitStructs[Nb_implicitfunction].zmin = (tmp.substr(0,position));
  553                 Nb_implicitfunction++;
  554                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  555                 ImplicitFct = tmp2;
  556             }
  557             else
  558             {
  559                 ImplicitStructs[Nb_implicitfunction].zmin  = ImplicitFct;
  560                 ImplicitFct ="";
  561             }
  562         }
  563         return Nb_implicitfunction;
  564     }
  565     else if(type == 6) //zmax
  566     {
  567         ImplicitStructs[0].zmax = ImplicitFct;
  568         while( ImplicitFct!= "")
  569         {
  570             if((position = ImplicitFct.find(";")) != std::string::npos)
  571             {
  572                 tmp = ImplicitFct;
  573                 ImplicitStructs[Nb_implicitfunction].zmax = (tmp.substr(0,position));
  574                 Nb_implicitfunction++;
  575                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  576                 ImplicitFct = tmp2;
  577             }
  578             else
  579             {
  580                 ImplicitStructs[Nb_implicitfunction].zmax  = ImplicitFct;
  581                 ImplicitFct ="";
  582             }
  583         }
  584         return Nb_implicitfunction;
  585     }
  586     else if(type == 8) //Cnd
  587     {
  588         ImplicitStructs[0].cnd = ImplicitFct;
  589         while( ImplicitFct!= "")
  590         {
  591             if((position = ImplicitFct.find(";")) != std::string::npos)
  592             {
  593                 tmp = ImplicitFct;
  594                 ImplicitStructs[Nb_implicitfunction].cnd = (tmp.substr(0,position));
  595                 Nb_implicitfunction++;
  596                 tmp2 = ImplicitFct.substr(position+1, ImplicitFct.length()-1);
  597                 ImplicitFct = tmp2;
  598             }
  599             else
  600             {
  601                 ImplicitStructs[Nb_implicitfunction].cnd  = ImplicitFct;
  602                 ImplicitFct ="";
  603             }
  604         }
  605         return Nb_implicitfunction;
  606     }
  607     return 0;
  608 }
  609 void Iso3D::ReinitVarTablesWhenMorphActiv(uint IsoIndex)
  610 {
  611     double vals[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  612     const uint limitX = masterthread->XYZgrid, limitY = masterthread->XYZgrid, limitZ = masterthread->XYZgrid;
  613     uint maxgridval = masterthread->GridVal;
  614     vals[3]         = masterthread->stepMorph;
  615     masterthread->xLocal2[IsoIndex*maxgridval]=(masterthread->x_Sup[IsoIndex]=masterthread->xSupParser[IsoIndex].Eval(vals));
  616     masterthread->yLocal2[IsoIndex*maxgridval]=(masterthread->y_Sup[IsoIndex]=masterthread->ySupParser[IsoIndex].Eval(vals));
  617     masterthread->zLocal2[IsoIndex*maxgridval]=(masterthread->z_Sup[IsoIndex]=masterthread->zSupParser[IsoIndex].Eval(vals));
  618     masterthread->x_Step[IsoIndex] = (masterthread->xLocal2[IsoIndex*maxgridval]-(masterthread->x_Inf[IsoIndex]=masterthread->xInfParser[IsoIndex].Eval(vals)))/(limitX-1);
  619     masterthread->y_Step[IsoIndex] = (masterthread->yLocal2[IsoIndex*maxgridval]-(masterthread->y_Inf[IsoIndex]=masterthread->yInfParser[IsoIndex].Eval(vals)))/(limitY-1);
  620     masterthread->z_Step[IsoIndex] = (masterthread->zLocal2[IsoIndex*maxgridval]-(masterthread->z_Inf[IsoIndex]=masterthread->zInfParser[IsoIndex].Eval(vals)))/(limitZ-1);
  621     for (uint i= 1; i < limitX; i++)
  622         masterthread->xLocal2[IsoIndex*maxgridval+i] = masterthread->xLocal2[IsoIndex*maxgridval+i-1] - masterthread->x_Step[IsoIndex];
  623     for (uint j= 1; j < limitY; j++)
  624         masterthread->yLocal2[IsoIndex*maxgridval+j] = masterthread->yLocal2[IsoIndex*maxgridval+j-1] - masterthread->y_Step[IsoIndex];
  625     for (uint k= 1; k < limitZ; k++)
  626         masterthread->zLocal2[IsoIndex*maxgridval+k] = masterthread->zLocal2[IsoIndex*maxgridval+k-1] - masterthread->z_Step[IsoIndex];
  627     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
  628     {
  629         for (uint k= 0; k < limitX; k++)
  630             workerthreads[nbthreads].xLocal2[IsoIndex*maxgridval+k] = masterthread->xLocal2[IsoIndex*maxgridval+k];
  631         for (uint k= 0; k < limitY; k++)
  632             workerthreads[nbthreads].yLocal2[IsoIndex*maxgridval+k] = masterthread->yLocal2[IsoIndex*maxgridval+k];
  633         for (uint k= 0; k < limitZ; k++)
  634             workerthreads[nbthreads].zLocal2[IsoIndex*maxgridval+k] = masterthread->zLocal2[IsoIndex*maxgridval+k];
  635     }
  636 }
  637 void IsoWorkerThread::VoxelEvaluation(uint IsoIndex)
  638 {
  639     double* vals;
  640     double* Res;
  641     uint maxgrscalemaxgr = GridVal*GridVal;
  642     const uint limitY = XYZgrid, limitZ = XYZgrid;
  643     uint I, J, IJK;
  644     uint id=0;
  645     uint nbX=OrignbX, nbY=OrignbY, nbZ=OrignbZ;
  646     uint nbstack=nbX*nbY*nbZ;
  647     uint Iindice=0, Jindice=0, Kindice=0, nbvar=8;
  648     int PreviousSignal=0;
  649     vals = new double[nbvar*nbstack];
  650     Res  = new double[nbstack];
  651     vals[3]    = stepMorph;
  652     //vals[4]    = double(MyIndex);
  653     uint taille=0;
  654     iStart = 0;
  655     iFinish = 0;
  656     for(uint i=0; i<XYZgrid; i++)
  657     {
  658         if((i% (WorkerThreadsNumber))  == MyIndex )
  659             taille += 1;
  660         if(MyIndex <= (i% (WorkerThreadsNumber)))
  661             iFinish  += 1;
  662     }
  663     iStart = iFinish - taille;
  664     for(uint l=0; l<nbstack; l++)
  665         vals[l*nbvar+3]= stepMorph;
  666     for(uint l=0; l<nbstack; l++)
  667         vals[l*nbvar+7]=GridVal;
  668     uint remX= (iFinish-iStart)%nbX;
  669     uint remY= limitY%nbY;
  670     uint remZ= limitZ%nbZ;
  671     uint Totalpoints=(iFinish-iStart)*limitY*limitZ;
  672     implicitFunctionParser[IsoIndex].AllocateStackMemory(Stack_Factor, nbvariables);
  673     for(uint i=iStart; i<iFinish; i+=nbX )
  674     {
  675         Iindice = i;
  676         nbY=OrignbY;
  677         if((remX>0) && ((Iindice+remX)==(iFinish)))
  678         {
  679             nbX = remX;
  680             i= iFinish;
  681         }
  682         I = Iindice*maxgrscalemaxgr;
  683         for(uint j=0; j<limitY; j+=nbY)
  684         {
  685             Jindice = j;
  686             nbZ=OrignbZ;
  687             if((remY>0) && ((Jindice+remY)==(limitY)))
  688             {
  689                 nbY = remY;
  690                 j= limitY;
  691             }
  692             J = I + Jindice*GridVal;
  693             for(uint k=0; k<limitZ; k+=nbZ)
  694             {
  695                 Kindice = k;
  696                 if((remZ>0) && ((Kindice+remZ)==(limitZ)))
  697                 {
  698                     nbZ = remZ;
  699                     k= limitZ;
  700                 }
  701                 nbstack = nbX*nbY*nbZ;
  702                 // X, Y and Z arrays construction:
  703                 for(uint l=0; l<nbstack; l++)
  704                 {
  705                     vals[l*nbvar  ]= xLocal2[IsoIndex*GridVal+Iindice+uint(l*nbX/nbstack)];
  706                     vals[l*nbvar+1]= yLocal2[IsoIndex*GridVal+Jindice+((uint(l/nbZ))%nbY)];
  707                     vals[l*nbvar+2]= zLocal2[IsoIndex*GridVal+Kindice+(l%nbZ)];
  708 
  709                     vals[l*nbvar+4]=Iindice+uint(l*nbX/nbstack);
  710                     vals[l*nbvar+5]=Jindice+((uint(l/nbZ))%nbY);
  711                     vals[l*nbvar+6]=Kindice+(l%nbZ);
  712                 }
  713                 IJK = J+Kindice;
  714                 double res = implicitFunctionParser[IsoIndex].Eval2(vals, nbvar, Res, nbstack);
  715                 if( abs(res - IF_FUNCT_ERROR) == 0.0)
  716                 {
  717                     for(uint l=0; l<nbstack; l++)
  718                         Res[l] = implicitFunctionParser[IsoIndex].Eval(&(vals[l*nbvar]));
  719                 }
  720                 else if( abs(res - DIVISION_BY_ZERO) == 0.0 || abs(res - VAR_OVERFLOW) == 0.0)
  721                 {
  722                     StopCalculations = true;
  723                 }
  724                 if(StopCalculations)
  725                     return;
  726                 uint p=0;
  727                 uint sect=0;
  728                 for(uint ii=0; ii<nbX; ii++)
  729                     for(uint jj=0; jj<nbY; jj++)
  730                         for(uint kk=0; kk<nbZ; kk++)
  731                         {
  732                             sect= IJK + (ii)*GridVal*GridVal+ (jj)*GridVal + (kk);
  733                             Results[sect] = Res[p];
  734                             GridVoxelVarPt[sect].Signature   = 0; // Signature initialisation
  735                             GridVoxelVarPt[sect].NbEdgePoint = 0; // No EdgePoint yet!
  736                             p++;
  737                         }
  738                 //Signal emission:
  739                 id+=nbstack;
  740                 if(MyIndex == 0 && activeMorph != 1)
  741                 {
  742                     signalVal = int((id*100)/Totalpoints);
  743                     if((signalVal - PreviousSignal) > 1 || id==Totalpoints)
  744                     {
  745                         PreviousSignal = signalVal;
  746                         emitMySignal();
  747                     }
  748                 }
  749             }
  750         }
  751     }
  752     delete[] vals;
  753     delete[] Res;
  754 }
  755 void Iso3D::ConstructIsoNormale(uint idx)
  756 {
  757     float val1, val2, val3, val4, val5, val6,
  758           pt1_x, pt1_y, pt1_z,
  759           pt2_x, pt2_y, pt2_z,
  760           pt3_x, pt3_y, pt3_z,
  761           scalar, n1, n2, n3;
  762     uint ThreeTimesI, IndexFirstPoint, IndexSecondPoint, IndexThirdPoint;
  763     NormOriginaltmpVector.resize(NormOriginaltmpVector.size()+ 3*NbTriangleIsoSurface);
  764     for(uint i = 0; i<NbTriangleIsoSurface; ++i)
  765     {
  766         ThreeTimesI      = i*3;
  767         IndexFirstPoint  = 10*IndexPolyTabVector[ThreeTimesI  +3*idx] + 7;
  768         IndexSecondPoint = 10*IndexPolyTabVector[ThreeTimesI+1+3*idx] + 7;
  769         IndexThirdPoint  = 10*IndexPolyTabVector[ThreeTimesI+2+3*idx] + 7;
  770         pt1_x= NormVertexTabVector[IndexFirstPoint  ];
  771         pt1_y= NormVertexTabVector[IndexFirstPoint+1];
  772         pt1_z= NormVertexTabVector[IndexFirstPoint+2];
  773         pt2_x= NormVertexTabVector[IndexSecondPoint  ];
  774         pt2_y= NormVertexTabVector[IndexSecondPoint+1];
  775         pt2_z= NormVertexTabVector[IndexSecondPoint+2];
  776         pt3_x= NormVertexTabVector[IndexThirdPoint   ];
  777         pt3_y= NormVertexTabVector[IndexThirdPoint+1 ];
  778         pt3_z= NormVertexTabVector[IndexThirdPoint+2 ];
  779         val1 = pt2_y - pt1_y;
  780         val2 = pt3_z - pt1_z;
  781         val3 = pt2_z - pt1_z;
  782         val4 = pt3_y - pt1_y;
  783         val5 = pt3_x - pt1_x;
  784         val6 = pt2_x - pt1_x;
  785         n1 = (val1*val2 - val3*val4);
  786         n2 = (val3*val5 - val6*val2);
  787         n3 = (val6*val4 - val1*val5);
  788         scalar = float(sqrt(n1*n1+n2*n2+n3*n3));
  789         if(scalar < float(0.0000001))  scalar  = float(0.0000001);
  790         (NormOriginaltmpVector[ThreeTimesI  ] = n1/scalar);
  791         (NormOriginaltmpVector[ThreeTimesI+1] = n2/scalar);
  792         (NormOriginaltmpVector[ThreeTimesI+2] = n3/scalar);
  793     }
  794 }
  795 void Iso3D::SaveIsoGLMap(uint indx)
  796 {
  797     uint IndexFirstPoint, IndexSecondPoint, IndexThirdPoint, ThreeTimesI;
  798     double scalar;
  799 /// Recalculate the normals so we have one for each Point (like Pov Mesh) :
  800     for (uint i=0; i < NbPointIsoMap ; i++)
  801     {
  802         ThreeTimesI = 10*i  + 4;
  803         NormVertexTabVector[ 10*NbVertexTmp +ThreeTimesI  ] = 0;
  804         NormVertexTabVector[ 10*NbVertexTmp +ThreeTimesI+1] = 0;
  805         NormVertexTabVector[ 10*NbVertexTmp +ThreeTimesI+2] = 0;
  806     }
  807     for(uint i = 0; i<NbTriangleIsoSurface; ++i)
  808     {
  809         ThreeTimesI   = i*3;
  810         IndexFirstPoint  = 10*IndexPolyTabVector[ThreeTimesI  +3*indx] + 4;
  811         IndexSecondPoint = 10*IndexPolyTabVector[ThreeTimesI+1+3*indx] + 4;
  812         IndexThirdPoint  = 10*IndexPolyTabVector[ThreeTimesI+2+3*indx] + 4;
  813         NormVertexTabVector[IndexFirstPoint  ] += NormOriginaltmpVector[ThreeTimesI  ];
  814         NormVertexTabVector[IndexFirstPoint+1] += NormOriginaltmpVector[ThreeTimesI+1];
  815         NormVertexTabVector[IndexFirstPoint+2] += NormOriginaltmpVector[ThreeTimesI+2];
  816         NormVertexTabVector[IndexSecondPoint  ] += NormOriginaltmpVector[ThreeTimesI  ];
  817         NormVertexTabVector[IndexSecondPoint+1] += NormOriginaltmpVector[ThreeTimesI+1];
  818         NormVertexTabVector[IndexSecondPoint+2] += NormOriginaltmpVector[ThreeTimesI+2];
  819         NormVertexTabVector[IndexThirdPoint  ]  += NormOriginaltmpVector[ThreeTimesI  ];
  820         NormVertexTabVector[IndexThirdPoint+1]  += NormOriginaltmpVector[ThreeTimesI+1];
  821         NormVertexTabVector[IndexThirdPoint+2]  += NormOriginaltmpVector[ThreeTimesI+2];
  822     }
  823 /// Normalisation
  824     uint idx;
  825     for (uint i=0; i < NbPointIsoMap  ; i++)
  826     {
  827         idx = 10*i + 4;
  828         scalar = double(sqrt((NormVertexTabVector[idx  ]*NormVertexTabVector[idx  ]) +
  829                              (NormVertexTabVector[idx+1]*NormVertexTabVector[idx+1]) +
  830                              (NormVertexTabVector[idx+2]*NormVertexTabVector[idx+2])));
  831         if(scalar < 0.000000001)  scalar = 0.000000001;
  832         NormVertexTabVector[idx  ] /= float(scalar);
  833         NormVertexTabVector[idx+1] /= float(scalar);
  834         NormVertexTabVector[idx+2] /= float(scalar);
  835     }
  836 }
  837 ErrorMessage IsoMasterThread::ParserIso()
  838 {
  839     double vals[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  840     initparser();
  841     //Evaluates defined constantes:
  842     if(constnotnull)
  843     {
  844         ConstSize = HowManyVariables(Const, 1);
  845         for(uint j=0; j<ConstSize; j++)
  846         {
  847             if ((stdError.iErrorIndex = Cstparser.Parse(Consts[j],"u")) >= 0)
  848             {
  849                 stdError.strError = Consts[j];
  850                 return stdError;
  851             }
  852             ConstValues.push_back(Cstparser.Eval(&vals[3]));
  853             Cstparser.AddConstant(ConstNames[j], ConstValues[j]);
  854         }
  855     }
  856     else
  857     {
  858         ConstSize = 0;
  859     }
  860     if(functnotnull)
  861     {
  862         FunctSize = HowManyVariables(Funct, 2);
  863         for(uint i=0; i<FunctSize; i++)
  864         {
  865             for(uint j=0; j<ConstSize; j++)
  866             {
  867                 Fct[i].AddConstant(ConstNames[j], ConstValues[j]);
  868             }
  869             //Add predefined constatnts:
  870             for(uint k=0; k<Nb_Sliders; k++)
  871             {
  872                 Fct[i].AddConstant(SliderNames[k], SliderValues[k]);
  873             }
  874         }
  875         for(uint i=0; i<FunctSize; i++)
  876         {
  877             for(uint j=0; j<i; j++)
  878                 if( (UsedFunct2[i*FunctSize+j]=(Functs[i].find(FunctNames[j]) != std::string::npos)))
  879                     Fct[i].AddFunction(FunctNames[j], Fct[j]);
  880             if ((stdError.iErrorIndex = Fct[i].Parse(Functs[i],"x,y,z,t"))>=0)
  881             {
  882                 stdError.strError = Functs[i];
  883                 return stdError;
  884             }
  885             Fct[i].AllocateStackMemory(Stack_Factor, nbvariables);
  886         }
  887     }
  888     else
  889     {
  890         FunctSize =0;
  891     }
  892     if(rgbtnotnull)
  893     {
  894         RgbtSize = HowManyVariables(Rgbt, 3);
  895         for(uint i=0; i<RgbtSize; i++)
  896         {
  897             for(uint j=0; j<ConstSize; j++)
  898             {
  899                 RgbtParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  900                 RgbtParser[i].AddConstant("Lacunarity", double(Lacunarity));
  901                 RgbtParser[i].AddConstant("Gain", double(Gain));
  902                 RgbtParser[i].AddConstant("Octaves", double(Octaves));
  903                 RgbtParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
  904                 RgbtParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
  905                 RgbtParser[i].AddFunction("MarbleP",MarblePerlin, 4);
  906             }
  907             //Add predefined constatnts:
  908             for(uint k=0; k<Nb_Sliders; k++)
  909             {
  910                 RgbtParser[i].AddConstant(SliderNames[k], SliderValues[k]);
  911             }
  912         }
  913     }
  914     else
  915     {
  916         RgbtSize =0;
  917     }
  918     //For Solid Texture :
  919     if(vrgbtnotnull)
  920     {
  921         VRgbtSize = HowManyVariables(VRgbt, 4);
  922         for(uint j=0; j<ConstSize; j++)
  923         {
  924             GradientParser->AddConstant(ConstNames[j], ConstValues[j]);
  925             //Add predefined constatnts:
  926             for(uint k=0; k<Nb_Sliders; k++)
  927                 GradientParser->AddConstant(SliderNames[k], SliderValues[k]);
  928         }
  929         GradientParser->AddConstant("Lacunarity", Lacunarity);
  930         GradientParser->AddConstant("Gain", Gain);
  931         GradientParser->AddConstant("Octaves", Octaves);
  932         GradientParser->AddFunction("NoiseW",TurbulenceWorley, 6);
  933         GradientParser->AddFunction("NoiseP",TurbulencePerlin, 6);
  934         GradientParser->AddFunction("MarbleP",MarblePerlin, 4);
  935         for(uint i=0; i<VRgbtSize; i++)
  936         {
  937             for(uint j=0; j<ConstSize; j++)
  938             {
  939                 VRgbtParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  940                 //Add predefined constatnts:
  941                 for(uint k=0; k<Nb_Sliders; k++)
  942                     VRgbtParser[i].AddConstant(SliderNames[k], SliderValues[k]);
  943                 VRgbtParser[i].AddConstant("Lacunarity", Lacunarity);
  944                 VRgbtParser[i].AddConstant("Gain", Gain);
  945                 VRgbtParser[i].AddConstant("Octaves", Octaves);
  946                 VRgbtParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
  947                 VRgbtParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
  948                 VRgbtParser[i].AddFunction("MarbleP",MarblePerlin, 4);
  949             }
  950         }
  951     }
  952     else
  953     {
  954         VRgbtSize =0;
  955     }
  956     if(Noise != "")
  957     {
  958         for(uint j=0; j<ConstSize; j++)
  959             NoiseParser->AddConstant(ConstNames[j], ConstValues[j]);
  960         NoiseParser->AddConstant("Lacunarity", Lacunarity);
  961         NoiseParser->AddConstant("Gain", Gain);
  962         NoiseParser->AddConstant("Octaves", Octaves);
  963         //Add predefined constatnts:
  964         for(uint k=0; k<Nb_Sliders; k++)
  965         {
  966             NoiseParser->AddConstant(SliderNames[k], SliderValues[k]);
  967         }
  968     }
  969     //ImplicitFunction is composed of more than one isosurface:
  970     HowManyIsosurface(ImplicitFunction, 0);
  971     HowManyIsosurface(XlimitInf, 1);
  972     HowManyIsosurface(XlimitSup, 2);
  973     HowManyIsosurface(YlimitInf, 3);
  974     HowManyIsosurface(YlimitSup, 4);
  975     HowManyIsosurface(ZlimitInf, 5);
  976     HowManyIsosurface(ZlimitSup, 6);
  977     HowManyIsosurface(Grid, 7);
  978     if(cndnotnull)
  979     {
  980         ParisoCondition = 1;
  981         HowManyIsosurface(Condition, 8);
  982     }
  983     else
  984         ParisoCondition = -1;
  985     //Add defined constantes:
  986     for(uint i=0; i<componentsNumber; i++)
  987     {
  988         for(uint j=0; j<ConstSize; j++)
  989         {
  990             implicitFunctionParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  991             if(cndnotnull)
  992                 ParisoConditionParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  993             xSupParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  994             ySupParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  995             zSupParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  996             xInfParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  997             yInfParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  998             zInfParser[i].AddConstant(ConstNames[j], ConstValues[j]);
  999         }
 1000         //Add predefined constatnts:
 1001         for(uint k=0; k<Nb_Sliders; k++)
 1002         {
 1003             implicitFunctionParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1004             if(cndnotnull)
 1005                 ParisoConditionParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1006             xSupParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1007             ySupParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1008             zSupParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1009             xInfParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1010             yInfParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1011             zInfParser[i].AddConstant(SliderNames[k], SliderValues[k]);
 1012         }
 1013     }
 1014     // Add defined functions :
 1015     if(rgbtnotnull)
 1016         for(int i=0; i<4; i++)
 1017             for(uint j=0; j<FunctSize; j++)
 1018             {
 1019                 RgbtParser[i].AddFunction(FunctNames[j], Fct[j]);
 1020                 RgbtParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
 1021                 RgbtParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
 1022                 RgbtParser[i].AddFunction("MarbleP",MarblePerlin, 4);
 1023             }
 1024     // Add defined functions :
 1025     if(vrgbtnotnull)
 1026     {
 1027         for(uint j=0; j<FunctSize; j++)
 1028         {
 1029             GradientParser->AddFunction(FunctNames[j], Fct[j]);
 1030             GradientParser->AddFunction("NoiseW",TurbulenceWorley, 6);
 1031             GradientParser->AddFunction("NoiseP",TurbulencePerlin, 6);
 1032             GradientParser->AddFunction("MarbleP",MarblePerlin, 4);
 1033         }
 1034 
 1035         for(int i=0; i<4; i++)
 1036             for(uint j=0; j<FunctSize; j++)
 1037             {
 1038                 VRgbtParser[i].AddFunction(FunctNames[j], Fct[j]);
 1039                 VRgbtParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
 1040                 VRgbtParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
 1041                 VRgbtParser[i].AddFunction("MarbleP",MarblePerlin, 4);
 1042             }
 1043     }
 1044     //delete NoiseFunction;
 1045     //NoiseFunction = new CellNoise();
 1046     for(uint i=0; i<componentsNumber; i++)
 1047     {
 1048         for(uint j=0; j<FunctSize; j++)
 1049         {
 1050             if((UsedFunct[i*FunctSize+j]=(ImplicitStructs[i].fxyz.find(FunctNames[j]) != std::string::npos)))
 1051             {
 1052                 implicitFunctionParser[i].AddFunction(FunctNames[j], Fct[j]);
 1053                 ParisoConditionParser[i].AddFunction(FunctNames[j], Fct[j]);
 1054             }
 1055         }
 1056         implicitFunctionParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
 1057         implicitFunctionParser[i].AddFunction("fhelix1",fhelix1, 10);
 1058         implicitFunctionParser[i].AddFunction("fhelix2",fhelix2, 10);
 1059         implicitFunctionParser[i].AddFunction("f_hex_y",f_hex_y, 4);
 1060         implicitFunctionParser[i].AddFunction("p_skeletal_int",p_skeletal_int, 3);
 1061         implicitFunctionParser[i].AddFunction("fmesh",fmesh, 8);
 1062         implicitFunctionParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
 1063         implicitFunctionParser[i].AddFunction("MarbleP",MarblePerlin, 4);
 1064     }
 1065     NoiseParser->AddFunction("NoiseW",TurbulenceWorley, 6);
 1066     NoiseParser->AddFunction("NoiseP",TurbulencePerlin, 6);
 1067     NoiseParser->AddFunction("MarbleP",MarblePerlin, 4);
 1068     return ParseExpression();
 1069 }
 1070 ErrorMessage IsoMasterThread::ParseExpression()
 1071 {
 1072     double vals[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 1073     uint limitX = XYZgrid, limitY = XYZgrid, limitZ = XYZgrid;
 1074 
 1075     if(AllComponentTraited)
 1076     {
 1077         stepMorph += pace;
 1078     }
 1079     vals[3] = stepMorph;
 1080     // Parse
 1081     if(rgbtnotnull && RgbtSize == 4)
 1082         for(uint i=0; i<RgbtSize; i++)
 1083             if ((stdError.iErrorIndex = RgbtParser[i].Parse(Rgbts[i],"x,y,z,t,cmpId,indx,x_step,y_step,z_step,max_ijk,x_sup,y_sup,z_sup,x_inf,y_inf,z_inf")) >= 0)
 1084             {
 1085                 stdError.strError = Rgbts[i];
 1086                 return stdError;
 1087             }
 1088     // Parse
 1089     if(vrgbtnotnull && (VRgbtSize % 5) ==0)
 1090     {
 1091         if ((stdError.iErrorIndex = GradientParser->Parse(Gradient,"x,y,z,t,cmpId,indx,x_step,y_step,z_step,max_ijk,x_sup,y_sup,z_sup,x_inf,y_inf,z_inf")) >= 0)
 1092         {
 1093             stdError.strError = Gradient;
 1094             return stdError;
 1095         }
 1096         for(uint i=0; i<VRgbtSize; i++)
 1097             if ((stdError.iErrorIndex = VRgbtParser[i].Parse(VRgbts[i],"x,y,z,t,cmpId,indx,x_step,y_step,z_step,max_ijk,x_sup,y_sup,z_sup,x_inf,y_inf,z_inf")) >= 0)
 1098             {
 1099                 stdError.strError = VRgbts[i];
 1100                 return stdError;
 1101             }
 1102     }
 1103     if(Noise != "")
 1104         if ((stdError.iErrorIndex = NoiseParser->Parse(Noise,"x,y,z,t,cmpId,indx,x_step,y_step,z_step,max_ijk,x_sup,y_sup,z_sup,x_inf,y_inf,z_inf")) >= 0)
 1105         {
 1106             stdError.strError = Noise;
 1107             return stdError;
 1108         }
 1109     for(uint i=0; i<componentsNumber; i++)
 1110     {
 1111         if ((stdError.iErrorIndex = implicitFunctionParser[i].Parse(ImplicitStructs[i].fxyz,"x,y,z,t,i_indx,j_indx,k_indx,max_ijk")) >= 0)
 1112         {
 1113             stdError.strError = ImplicitStructs[i].fxyz;
 1114             return stdError;
 1115         }
 1116         if(cndnotnull && (ImplicitStructs[i].cnd!=""))
 1117         {
 1118             if ((stdError.iErrorIndex = ParisoConditionParser[i].Parse(ImplicitStructs[i].cnd,"x,y,z,t")) >= 0)
 1119             {
 1120                 stdError.strError = ImplicitStructs[i].cnd;
 1121                 return stdError;
 1122             }
 1123         }
 1124         if ((stdError.iErrorIndex = xSupParser[i].Parse(ImplicitStructs[i].xmax, "x,y,z,t")) >= 0)
 1125         {
 1126             stdError.strError = ImplicitStructs[i].xmax;
 1127             return stdError;
 1128         }
 1129         if ((stdError.iErrorIndex = ySupParser[i].Parse(ImplicitStructs[i].ymax, "x,y,z,t")) >= 0)
 1130         {
 1131             stdError.strError = ImplicitStructs[i].ymax;
 1132             return stdError;
 1133         }
 1134         if ((stdError.iErrorIndex = zSupParser[i].Parse(ImplicitStructs[i].zmax, "x,y,z,t")) >= 0)
 1135         {
 1136             stdError.strError = ImplicitStructs[i].zmax;
 1137             return stdError;
 1138         }
 1139         if ((stdError.iErrorIndex = xInfParser[i].Parse(ImplicitStructs[i].xmin, "x,y,z,t")) >= 0)
 1140         {
 1141             stdError.strError = ImplicitStructs[i].xmin;
 1142             return stdError;
 1143         }
 1144         if ((stdError.iErrorIndex = yInfParser[i].Parse(ImplicitStructs[i].ymin, "x,y,z,t")) >= 0)
 1145         {
 1146             stdError.strError = ImplicitStructs[i].ymin;
 1147             return stdError;
 1148         }
 1149         if ((stdError.iErrorIndex = zInfParser[i].Parse(ImplicitStructs[i].zmin, "x,y,z,t")) >= 0)
 1150         {
 1151             stdError.strError = ImplicitStructs[i].zmin;
 1152             return stdError;
 1153         }
 1154     }
 1155     for(uint IsoIndex=0; IsoIndex<componentsNumber; IsoIndex++)
 1156     {
 1157         if(gridnotnull)
 1158         {
 1159             limitX = limitY = limitZ = grid[IsoIndex];
 1160         }
 1161         xLocal2[IsoIndex*GridVal]=(x_Sup[IsoIndex]=xSupParser[IsoIndex].Eval(vals));
 1162         yLocal2[IsoIndex*GridVal]=(y_Sup[IsoIndex]=ySupParser[IsoIndex].Eval(vals));
 1163         zLocal2[IsoIndex*GridVal]=(z_Sup[IsoIndex]=zSupParser[IsoIndex].Eval(vals));
 1164         x_Step[IsoIndex] = (xLocal2[IsoIndex*GridVal] - (x_Inf[IsoIndex]=xInfParser[IsoIndex].Eval(vals)))/(limitX-1);
 1165         y_Step[IsoIndex] = (yLocal2[IsoIndex*GridVal] - (y_Inf[IsoIndex]=yInfParser[IsoIndex].Eval(vals)))/(limitY-1);
 1166         z_Step[IsoIndex] = (zLocal2[IsoIndex*GridVal] - (z_Inf[IsoIndex]=zInfParser[IsoIndex].Eval(vals)))/(limitZ-1);
 1167         for (uint i= 1; i < limitX; i++) xLocal2[IsoIndex*GridVal+i] = xLocal2[IsoIndex*GridVal+i-1] - x_Step[IsoIndex];
 1168         for (uint j= 1; j < limitY; j++) yLocal2[IsoIndex*GridVal+j] = yLocal2[IsoIndex*GridVal+j-1] - y_Step[IsoIndex];
 1169         for (uint k= 1; k < limitZ; k++) zLocal2[IsoIndex*GridVal+k] = zLocal2[IsoIndex*GridVal+k-1] - z_Step[IsoIndex];
 1170     }
 1171     return stdError;
 1172 }
 1173 void IsoMasterThread::DeleteMasterParsers()
 1174 {
 1175     if(ParsersAllocated)
 1176     {
 1177         delete[] implicitFunctionParser;
 1178         delete[] xSupParser;
 1179         delete[] ySupParser;
 1180         delete[] zSupParser;
 1181         delete[] xInfParser;
 1182         delete[] yInfParser;
 1183         delete[] zInfParser;
 1184         delete[] Fct;
 1185         delete[] UsedFunct;
 1186         delete[] UsedFunct2;
 1187         delete[] ParisoConditionParser;
 1188         delete[] VRgbtParser;
 1189         delete[] RgbtParser;
 1190         delete GradientParser;
 1191         delete NoiseParser;
 1192         ParsersAllocated = false;
 1193     }
 1194     ImplicitStructs.clear();
 1195     xLocal2.clear();
 1196     yLocal2.clear();
 1197     zLocal2.clear();
 1198     x_Step.clear();
 1199     y_Step.clear();
 1200     z_Step.clear();
 1201     x_Sup.clear();
 1202     y_Sup.clear();
 1203     z_Sup.clear();
 1204     x_Inf.clear();
 1205     y_Inf.clear();
 1206     z_Inf.clear();
 1207     SliderValues.clear();
 1208     SliderNames.clear();
 1209     Consts.clear();
 1210     ConstNames.clear();
 1211     ConstValues.clear();
 1212     Rgbts.clear();
 1213     RgbtNames.clear();
 1214     VRgbts.clear();
 1215     VRgbtNames.clear();
 1216     Functs.clear();
 1217     FunctNames.clear();
 1218 }
 1219 void IsoWorkerThread::DeleteWorkerParsers()
 1220 {
 1221     if(ParsersAllocated)
 1222     {
 1223         delete[] implicitFunctionParser;
 1224         delete[] Fct;
 1225         ParsersAllocated = false;
 1226     }
 1227 }
 1228 void IsoMasterThread::InitMasterParsers()
 1229 {
 1230     for(uint i=0; i<componentsNumber; i++)
 1231     {
 1232         implicitFunctionParser[i].AddConstant("pi", PI);
 1233         implicitFunctionParser[i].AddConstant("ThreadId", MyIndex);
 1234         ParisoConditionParser[i].AddConstant("pi", PI);
 1235         xSupParser[i].AddConstant("pi", PI);
 1236         ySupParser[i].AddConstant("pi", PI);
 1237         zSupParser[i].AddConstant("pi", PI);
 1238         xInfParser[i].AddConstant("pi", PI);
 1239         yInfParser[i].AddConstant("pi", PI);
 1240         zInfParser[i].AddConstant("pi", PI);
 1241     }
 1242     NoiseParser->AddConstant("pi", PI);
 1243     NoiseParser->AddConstant("Lacunarity", Lacunarity);
 1244     NoiseParser->AddConstant("Gain", Gain);
 1245     NoiseParser->AddConstant("Octaves", Octaves);
 1246     for(uint i=0; i<FunctSize; i++)
 1247     {
 1248         Fct[i].AddConstant("pi", PI);
 1249         Fct[i].AddConstant("ThreadId", MyIndex);
 1250         Fct[i].AddFunction("CmpId",CurrentIsoCmpId, 1);
 1251         Fct[i].AddFunction("NoiseW",TurbulenceWorley, 6);
 1252         Fct[i].AddFunction("fhelix1",fhelix1, 10);
 1253         Fct[i].AddFunction("fhelix2",fhelix2, 10);
 1254         Fct[i].AddFunction("f_hex_y",f_hex_y, 4);
 1255         Fct[i].AddFunction("p_skeletal_int",p_skeletal_int, 3);
 1256         Fct[i].AddFunction("fmesh",fmesh, 8);
 1257         Fct[i].AddFunction("NoiseP",TurbulencePerlin, 6);
 1258         Fct[i].AddFunction("MarbleP",MarblePerlin, 4);
 1259     }
 1260     for(uint i=0; i< RgbtSize; i++)
 1261     {
 1262         RgbtParser[i].AddConstant("pi", PI);
 1263         RgbtParser[i].AddConstant("Lacunarity", Lacunarity);
 1264         RgbtParser[i].AddConstant("Gain", Gain);
 1265         RgbtParser[i].AddConstant("Octaves", Octaves);
 1266         RgbtParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
 1267         RgbtParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
 1268         RgbtParser[i].AddFunction("MarbleP",MarblePerlin, 4);
 1269     }
 1270     for(uint i=0; i<VRgbtSize; i++)
 1271     {
 1272         VRgbtParser[i].AddConstant("pi", PI);
 1273         VRgbtParser[i].AddConstant("Lacunarity", Lacunarity);
 1274         VRgbtParser[i].AddConstant("Gain", Gain);
 1275         VRgbtParser[i].AddConstant("Octaves", Octaves);
 1276         VRgbtParser[i].AddFunction("NoiseW",TurbulenceWorley, 6);
 1277         VRgbtParser[i].AddFunction("NoiseP",TurbulencePerlin, 6);
 1278         VRgbtParser[i].AddFunction("MarbleP",MarblePerlin, 4);
 1279     }
 1280     GradientParser->AddConstant("pi", PI);
 1281     GradientParser->AddConstant("Lacunarity", Lacunarity);
 1282     GradientParser->AddConstant("Gain", Gain);
 1283     GradientParser->AddConstant("Octaves", Octaves);
 1284     GradientParser->AddFunction("NoiseW",TurbulenceWorley, 6);
 1285     GradientParser->AddFunction("NoiseP",TurbulencePerlin, 6);
 1286     GradientParser->AddFunction("MarbleP",MarblePerlin, 4);
 1287 }
 1288 void IsoMasterThread::AllocateMasterParsers()
 1289 {
 1290     if(!ParsersAllocated)
 1291     {
 1292         implicitFunctionParser = new FunctionParser[componentsNumber];
 1293         xSupParser = new FunctionParser[componentsNumber];
 1294         ySupParser = new FunctionParser[componentsNumber];
 1295         zSupParser = new FunctionParser[componentsNumber];
 1296         xInfParser = new FunctionParser[componentsNumber];
 1297         yInfParser = new FunctionParser[componentsNumber];
 1298         zInfParser = new FunctionParser[componentsNumber];
 1299         ParisoConditionParser = new FunctionParser[componentsNumber];
 1300         ImplicitStructs.resize(componentsNumber);
 1301         xLocal2.resize(GridVal*componentsNumber);
 1302         yLocal2.resize(GridVal*componentsNumber);
 1303         zLocal2.resize(GridVal*componentsNumber);
 1304         x_Step.resize(componentsNumber);
 1305         y_Step.resize(componentsNumber);
 1306         z_Step.resize(componentsNumber);
 1307         x_Sup.resize(componentsNumber);
 1308         y_Sup.resize(componentsNumber);
 1309         z_Sup.resize(componentsNumber);
 1310         x_Inf.resize(componentsNumber);
 1311         y_Inf.resize(componentsNumber);
 1312         z_Inf.resize(componentsNumber);
 1313         if(!functnotnull)
 1314             FunctSize = 0;
 1315         Fct          = new FunctionParser[FunctSize];
 1316         UsedFunct    = new bool[4*componentsNumber*FunctSize];
 1317         UsedFunct2   = new bool[FunctSize*FunctSize];
 1318         vectnotnull? nbvariables=vect[0] : nbvariables=0;
 1319         rgbtnotnull ?
 1320         RgbtParser = new FunctionParser[(RgbtSize = 4)] :
 1321         RgbtParser = new FunctionParser[(RgbtSize = 0)];
 1322         vrgbtnotnull ?
 1323         VRgbtParser = new FunctionParser[VRgbtSize] :
 1324         VRgbtParser = new FunctionParser[(VRgbtSize = 0)];
 1325         if(constnotnull)
 1326             ConstSize=0;
 1327         GradientParser = new FunctionParser;
 1328         NoiseParser = new FunctionParser;
 1329         ParsersAllocated = true;
 1330     }
 1331 }
 1332 void IsoWorkerThread::AllocateParsersForWorkerThread(uint nbcomp, uint nbfunct)
 1333 {
 1334     if(!ParsersAllocated)
 1335     {
 1336         implicitFunctionParser = new FunctionParser[nbcomp];
 1337         Fct = new FunctionParser[nbfunct];
 1338         ParsersAllocated = true;
 1339     }
 1340 }
 1341 void IsoMasterThread::initgrid()
 1342 {
 1343     GridVal = XYZgrid;
 1344     if(gridnotnull)
 1345         for(uint fctnb= 0; fctnb< componentsNumber; fctnb++)
 1346             GridVal = std::max(GridVal, grid[fctnb]);
 1347 }
 1348 void IsoMasterThread::initparser()
 1349 {
 1350     DeleteMasterParsers();
 1351     initgrid();
 1352     AllocateMasterParsers();
 1353     InitMasterParsers();
 1354 }
 1355 void IsoWorkerThread::IsoCompute(uint fctnb)
 1356 {
 1357     (fctnb == 0) ? AllComponentTraited = true : AllComponentTraited = false;
 1358     VoxelEvaluation(fctnb);
 1359 }
 1360 void Iso3D::stopcalculations(bool calculation)
 1361 {
 1362     StopCalculations = calculation;
 1363     masterthread->StopCalculations = calculation;
 1364     for(uint nbthreads=0; nbthreads+1< WorkerThreadsNumber; nbthreads++)
 1365         workerthreads[nbthreads].StopCalculations = StopCalculations;
 1366 }
 1367 void Iso3D::copycomponent(struct ComponentInfos* copy, struct ComponentInfos* origin)
 1368 {
 1369     copy->ParisoTriangle = origin->ParisoTriangle;
 1370     copy->ParisoVertex          = origin->ParisoVertex;
 1371     copy->NbComponentsType    = origin->NbComponentsType;
 1372     copy->ParisoCurrentComponentIndex = origin->ParisoCurrentComponentIndex;
 1373     copy->ParisoNbComponents          = origin->ParisoNbComponents;
 1374     copy->Interleave                  = origin->Interleave;
 1375     copy->pariso                      = origin->pariso;
 1376     copy->updateviewer                = origin->updateviewer;
 1377     copy->ThereisCND                  = origin->ThereisCND;
 1378     copy->ParisoCondition             = origin->ParisoCondition;
 1379     copy->ThereisRGBA                 = origin->ThereisRGBA;
 1380     copy->NbTrianglesVerifyCND        = origin->NbTrianglesVerifyCND;
 1381     copy->NbTrianglesNoCND            = origin->NbTrianglesNoCND;
 1382     copy->NbTrianglesNotVerifyCND     = origin->NbTrianglesNotVerifyCND;
 1383     copy->NbTrianglesBorderCND        = origin->NbTrianglesBorderCND;
 1384     for(int i=0; i<2; i++)
 1385     {
 1386         copy->NoiseParam[i]  = origin->NoiseParam[i];
 1387     }
 1388 }
 1389 void Iso3D::IsoBuild (
 1390     float **NormVertexTabVectorPt,
 1391     uint **IndexPolyTabPt,
 1392     unsigned   int *PolyNumber,
 1393     uint *VertexNumberpt,
 1394     uint **IndexPolyTabMinPt,
 1395     unsigned  int *VertxNumber,
 1396     unsigned  int *MinimPolySize,
 1397     struct ComponentInfos * componentsPt
 1398 )
 1399 {
 1400     uint NbTriangleIsoSurfaceTmp, PreviousGridVal=masterthread->XYZgrid;
 1401     NbPointIsoMap= 0;
 1402     NbVertexTmp = NbTriangleIsoSurfaceTmp = 0;
 1403     componentsPt->updateviewer= false;
 1404     clear(components);
 1405     components->ParisoCondition = componentsPt->ParisoCondition;
 1406     if(componentsPt->pariso && componentsPt->ParisoCurrentComponentIndex>0)
 1407     {
 1408         NbVertexTmp = uint(NormVertexTabVector.size()/10);
 1409         NbTriangleIsoSurfaceTmp = uint(IndexPolyTabVector.size()/3);
 1410         copycomponent(components, componentsPt);
 1411     }
 1412     else
 1413     {
 1414         if(componentsPt->pariso)
 1415         {
 1416             components->pariso = true;
 1417             components->ParisoCurrentComponentIndex = componentsPt->ParisoCurrentComponentIndex;
 1418             components->ParisoNbComponents = componentsPt->ParisoNbComponents;
 1419         }
 1420         NormVertexTabVector.clear();
 1421         NormVertexTabVector.shrink_to_fit();
 1422         IndexPolyTabVector.clear();
 1423         IndexPolyTabVector.shrink_to_fit();
 1424         IndexPolyTabMinVector.clear();
 1425         IndexPolyTabMinVector.shrink_to_fit();
 1426     }
 1427     NormOriginaltmpVector.clear();
 1428     NormOriginaltmpVector.shrink_to_fit();
 1429     //*****//
 1430     uint maxx = std::max(masterthread->XYZgrid, masterthread->GridVal);
 1431     if(masterthread->gridnotnull)
 1432         for(uint fctnb= 0; fctnb< masterthread->componentsNumber; fctnb++)
 1433             maxx = std::max(maxx, masterthread->grid[fctnb]);
 1434     masterthread->GridVal = maxx;
 1435     for(uint nbthreads=0; nbthreads+1<WorkerThreadsNumber; nbthreads++)
 1436     {
 1437         workerthreads[nbthreads].GridVal = masterthread->GridVal;
 1438     }
 1439     if(GridVoxelVarPt != nullptr)
 1440         delete[] GridVoxelVarPt;
 1441     if(Results != nullptr)
 1442         delete[] Results;
 1443     try
 1444       {
 1445         GridVoxelVarPt = new Voxel[masterthread->GridVal*masterthread->GridVal*masterthread->GridVal];
 1446       }
 1447       catch (std::bad_alloc& e)
 1448       {
 1449         messageerror = MEM_OVERFLOW;
 1450         emitErrorSignal();
 1451         return;
 1452       }
 1453     Results = new (std::nothrow) double[masterthread->GridVal*masterthread->GridVal*masterthread->GridVal];
 1454     if (!Results)
 1455     {
 1456         messageerror = MEM_OVERFLOW;
 1457         emitErrorSignal();
 1458         return;
 1459     }
 1460     components->NbComponentsType.push_back(masterthread->componentsNumber);
 1461     stopcalculations(false);
 1462     if(masterthread->activeMorph != 1)
 1463     {
 1464         times.restart();
 1465     }
 1466     // generate Isosurface for all the implicit formulas
 1467     for(uint fctnb= 0; fctnb< masterthread->componentsNumber; fctnb++)
 1468     {
 1469         if(masterthread->activeMorph != 1)
 1470         {
 1471             message = QString("1) Cmp:"+QString::number(fctnb+1)+"/"+QString::number(masterthread->componentsNumber)+"==> Math calculation");
 1472             emitUpdateMessageSignal();
 1473         }
 1474         if(masterthread->gridnotnull)
 1475             Setgrid(masterthread->grid[fctnb]);
 1476         IsoComponentId = fctnb;
 1477         masterthread->CurrentComponent = fctnb;
 1478         for(uint nbthreads=0; nbthreads+1 < WorkerThreadsNumber; nbthreads++)
 1479             workerthreads[nbthreads].CurrentComponent = fctnb;
 1480         for(uint nbthreads=0; nbthreads+1 < WorkerThreadsNumber; nbthreads++)
 1481             workerthreads[nbthreads].stepMorph = masterthread->stepMorph;
 1482         if(masterthread->activeMorph == 1)
 1483         {
 1484             if(fctnb == 0)
 1485             {
 1486                 //This code is to ensure that stepmorph values are the same accross all threads
 1487                 masterthread->stepMorph += masterthread->pace;
 1488                 for(uint nbthreads=0; nbthreads+1 < WorkerThreadsNumber; nbthreads++)
 1489                     workerthreads[nbthreads].stepMorph = masterthread->stepMorph;
 1490             }
 1491             // Recalculate some tables values:
 1492             ReinitVarTablesWhenMorphActiv(fctnb);
 1493         }
 1494         masterthread->start();
 1495         for(uint nbthreads=0; nbthreads+1 < WorkerThreadsNumber; nbthreads++)
 1496             workerthreads[nbthreads].start();
 1497         masterthread->wait();
 1498         for(uint nbthreads=0; nbthreads+1 < WorkerThreadsNumber; nbthreads++)
 1499             workerthreads[nbthreads].wait();
 1500         bool Stop = masterthread->StopCalculations;
 1501         for(uint nbthreads=0; nbthreads+1 < WorkerThreadsNumber; nbthreads++)
 1502             Stop = Stop || workerthreads[nbthreads].StopCalculations;
 1503         if(StopCalculations || Stop)
 1504         {
 1505             Setgrid(PreviousGridVal);
 1506             delete[] GridVoxelVarPt;
 1507             GridVoxelVarPt = nullptr;
 1508             delete[] Results;
 1509             Results = nullptr;
 1510             return;
 1511         }
 1512         if(masterthread->activeMorph != 1)
 1513         {
 1514             message += QString(" ==> Mesh generation");
 1515             emitUpdateMessageSignal();
 1516         }
 1517         uint result = PointEdgeComputation(fctnb);
 1518         if(result == 0)
 1519         {
 1520             messageerror = VERTEX_TAB_MEM_OVERFLOW;
 1521             emitErrorSignal();
 1522             return;
 1523         }
 1524         SignatureComputation();
 1525         result = ConstructIsoSurface();
 1526         if(result == 0)
 1527         {
 1528             messageerror = TRIANGLES_TAB_MEM_OVERFLOW;
 1529             emitErrorSignal();
 1530             return;
 1531         }
 1532         ConstructIsoNormale(NbTriangleIsoSurfaceTmp);
 1533         SaveIsoGLMap(NbTriangleIsoSurfaceTmp);
 1534         NormOriginaltmpVector.clear();
 1535         result = SetMiniMmeshStruct();
 1536         if(result == 0)
 1537         {
 1538             messageerror = MINPOLY_TAB_MEM_OVERFLOW;
 1539             emitErrorSignal();
 1540             return;
 1541         }
 1542         // Save the Index:
 1543         components->ParisoTriangle.push_back(3*NbTriangleIsoSurfaceTmp); //save the starting position of this component
 1544         components->ParisoTriangle.push_back(NbTriangleIsoSurface);      //save the number of triangles of this component
 1545         components->ParisoVertex.push_back(NbVertexTmp);
 1546         components->ParisoVertex.push_back(NbVertexTmp + NbPointIsoMap -1);
 1547         // Save Number of Polys and vertex :
 1548         NbVertexTmp               += NbPointIsoMap;
 1549         NbTriangleIsoSurfaceTmp   += NbTriangleIsoSurface;
 1550     }
 1551     delete[] GridVoxelVarPt;
 1552     GridVoxelVarPt = nullptr;
 1553     delete[] Results;
 1554     Results = nullptr;
 1555     if(masterthread->activeMorph != 1)
 1556     {
 1557         message = QString("2) Mesh Processing");
 1558         emitUpdateMessageSignal();
 1559     }
 1560     //CND calculation for the triangles results:
 1561     uint result = CNDCalculation(NbTriangleIsoSurfaceTmp, components);
 1562     if(result == 0)
 1563     {
 1564         messageerror = CND_TAB_MEM_OVERFLOW;
 1565         emitErrorSignal();
 1566         return;
 1567     }
 1568     else if(result == 2)
 1569     {
 1570         messageerror = CND_POL_MEM_OVERFLOW;
 1571         emitErrorSignal();
 1572         return;
 1573     }
 1574     // Pigment, Texture and Noise :
 1575     if(masterthread->vrgbtnotnull)
 1576     {
 1577         components->ThereisRGBA.push_back(true);
 1578         components->NoiseParam[components->ParisoCurrentComponentIndex].NoiseType = 0; //Pigments
 1579         components->NoiseParam[components->ParisoCurrentComponentIndex].VRgbtParser = masterthread->VRgbtParser;
 1580         components->NoiseParam[components->ParisoCurrentComponentIndex].GradientParser = masterthread->GradientParser;
 1581         components->NoiseParam[components->ParisoCurrentComponentIndex].NoiseParser =  masterthread->NoiseParser;
 1582         components->NoiseParam[components->ParisoCurrentComponentIndex].Nb_vrgbts = masterthread->VRgbtSize;
 1583     }
 1584     else if(masterthread->RgbtSize >= 4)
 1585     {
 1586         components->ThereisRGBA.push_back(true);
 1587         components->NoiseParam[components->ParisoCurrentComponentIndex].NoiseType = 1; //Texture
 1588         components->NoiseParam[components->ParisoCurrentComponentIndex].RgbtParser = masterthread->RgbtParser;
 1589         components->NoiseParam[components->ParisoCurrentComponentIndex].NoiseParser =  masterthread->NoiseParser;
 1590     }
 1591     else
 1592     {
 1593         components->ThereisRGBA.push_back(false);
 1594         components->NoiseParam[components->ParisoCurrentComponentIndex].NoiseType = -1; //No Pigments or texture
 1595     }
 1596     if(masterthread->Noise == "")
 1597         components->NoiseParam[components->ParisoCurrentComponentIndex].NoiseShape = 0;
 1598     else
 1599         components->NoiseParam[components->ParisoCurrentComponentIndex].NoiseShape = 1;
 1600     CalculateColorsPoints(components, components->ThereisRGBA.size()-1);
 1601     *PolyNumber      = uint(IndexPolyTabVector.size());
 1602     *VertexNumberpt  = uint(NormVertexTabVector.size()/10);
 1603     *VertxNumber     = uint(IndexPolyTabMinVector.size());
 1604     if(masterthread->activeMorph != 1)
 1605     {
 1606         message = QString("Thr:"+QString::number(WorkerThreadsNumber)+"; Cmp:"+QString::number(masterthread->componentsNumber)+"; T="+QString::number(times.elapsed()/1000.0)+"s");
 1607         emitUpdateMessageSignal();
 1608     }
 1609 
 1610 
 1611     NormVertexTabVector.resize(NormVertexTabVector.size()+ (12+60+36)*10); // To add memory space to store the cube 12 vertices (three quads)
 1612 
 1613     uint startpl = 0;
 1614     uint actualpointindice=0;
 1615     for (uint i = 0; i < *VertxNumber; i++)
 1616     {
 1617         uint polysize = IndexPolyTabMinVector[startpl++];
 1618         for (uint j = 0; j < polysize; j++)
 1619         {
 1620             actualpointindice = IndexPolyTabMinVector[startpl];
 1621             IndexPolyTabVector.push_back(actualpointindice);
 1622             startpl++;
 1623         }
 1624         i += polysize;
 1625     }
 1626 
 1627 
 1628 
 1629 
 1630     IndexPolyTabMinVector2.clear();
 1631     for (uint i = 0; i < IndexPolyTabMinVector.size(); i++)
 1632     {
 1633         uint polysize = IndexPolyTabMinVector[i];
 1634         IndexPolyTabMinVector2.push_back(polysize);
 1635         i += polysize;
 1636     }
 1637 
 1638     *MinimPolySize = IndexPolyTabVector.size() - *PolyNumber;
 1639     *VertxNumber     = uint(IndexPolyTabMinVector2.size());
 1640     *IndexPolyTabMinPt = IndexPolyTabMinVector2.data();
 1641 
 1642     *NormVertexTabVectorPt = NormVertexTabVector.data();
 1643     *IndexPolyTabPt = IndexPolyTabVector.data();
 1644     copycomponent(componentsPt, components);
 1645     componentsPt->Interleave = true;
 1646     if(componentsPt->ParisoCurrentComponentIndex != (componentsPt->ParisoNbComponents-1))
 1647         componentsPt->ParisoCurrentComponentIndex += 1;
 1648     else
 1649         componentsPt->ParisoCurrentComponentIndex = 0;
 1650 
 1651     InitShowComponent(componentsPt);
 1652     Setgrid(PreviousGridVal);
 1653     componentsPt->updateviewer = true;
 1654 }
 1655 
 1656 void Iso3D::InitShowComponent(struct ComponentInfos *cpInfos)
 1657 {
 1658     cpInfos->ShowParIsoCmp.clear();
 1659     uint idx =0;
 1660     for(uint i=0; i<cpInfos->NbComponentsType.size(); i++)
 1661         idx+=cpInfos->NbComponentsType[i];
 1662     for(uint i=0; i<idx; i++)
 1663         cpInfos->ShowParIsoCmp.push_back(true);
 1664 }
 1665 
 1666 void Iso3D::Setgrid(uint NewGridVal)
 1667 {
 1668     if(masterthread->gridnotnull)
 1669     {
 1670         masterthread->XYZgrid  = NewGridVal;
 1671         for(uint th=0; th+1 < WorkerThreadsNumber; th++)
 1672             workerthreads[th].XYZgrid = NewGridVal;
 1673     }
 1674 }
 1675 void Iso3D::CalculateColorsPoints(struct ComponentInfos* comp, uint index)
 1676 {
 1677     uint cmpId=0, K=0;
 1678     double tmp,
 1679            *ValCol,
 1680            val[16];
 1681     ValCol = new double[masterthread->VRgbtSize];
 1682     val[3] = masterthread->stepMorph;
 1683     val[0] = val[1] = val[2] = 0.0;
 1684     if(comp->ThereisRGBA[index] == true &&  comp->NoiseParam[comp->ParisoCurrentComponentIndex].NoiseType == 0)
 1685     {
 1686         uint idx=0;
 1687         for(uint i=0; i < comp->NbComponentsType.size()-1; i++)
 1688             idx+=comp->NbComponentsType[i];
 1689         for(uint i= comp->ParisoVertex[2*idx]; i < NbVertexTmp; i++)
 1690         {
 1691             if((i >= uint(comp->ParisoVertex[2*(cmpId+idx)])))
 1692             {
 1693                 K = cmpId;
 1694                 if((masterthread->componentsNumber-1)>cmpId)
 1695                 {
 1696                     cmpId++;
 1697                 }
 1698             }
 1699             val[0]= double(NormVertexTabVector[i*10  + 7 ]);
 1700             val[1]= double(NormVertexTabVector[i*10  + 8 ]);
 1701             val[2]= double(NormVertexTabVector[i*10  + 9 ]);
 1702             val[4]= double(K);
 1703             if(masterthread->gridnotnull)
 1704             {
 1705                 val[5] = double(i);
 1706                 val[6] = masterthread->x_Step[K];
 1707                 val[7] = masterthread->y_Step[K];
 1708                 val[8] = masterthread->z_Step[K];
 1709                 val[9] = double(masterthread->grid[K]);
 1710                 val[10] = masterthread->x_Sup[K];
 1711                 val[11] = masterthread->y_Sup[K];
 1712                 val[12] = masterthread->z_Sup[K];
 1713                 val[13] = masterthread->x_Inf[K];
 1714                 val[14] = masterthread->y_Inf[K];
 1715                 val[15] = masterthread->z_Inf[K];
 1716             }
 1717             else
 1718             {
 1719                 val[5] = double(i);
 1720                 val[6] = masterthread->x_Step[0];
 1721                 val[7] = masterthread->y_Step[0];
 1722                 val[8] = masterthread->z_Step[0];
 1723                 val[9] = double(masterthread->GridVal);
 1724                 val[10] = masterthread->x_Sup[0];
 1725                 val[11] = masterthread->y_Sup[0];
 1726                 val[12] = masterthread->z_Sup[0];
 1727                 val[13] = masterthread->x_Inf[0];
 1728                 val[14] = masterthread->y_Inf[0];
 1729                 val[15] = masterthread->z_Inf[0];
 1730             }
 1731             for(uint i=0; i<masterthread->VRgbtSize; i++)
 1732             {
 1733                 ValCol[i] = masterthread->VRgbtParser[i].Eval(val);
 1734             }
 1735             if(masterthread->Noise != "")
 1736                 tmp  = masterthread->NoiseParser->Eval(val);
 1737             else
 1738                 tmp =1.0;
 1739             val[0]= tmp*double(NormVertexTabVector[i*10  + 7 ]);
 1740             val[1]= tmp*double(NormVertexTabVector[i*10  + 8 ]);
 1741             val[2]= tmp*double(NormVertexTabVector[i*10  + 9 ]);
 1742             tmp  = masterthread->GradientParser->Eval(val);
 1743             for (uint j=0; j < masterthread->VRgbtSize; j+=5)
 1744                 if(tmp < ValCol[j])
 1745                 {
 1746                     float fraction=0;
 1747                     if(j>=5 && (ValCol[j] != ValCol[j-5]))
 1748                     {
 1749                         fraction = (tmp-ValCol[j-5])/(ValCol[j]-ValCol[j-5]);
 1750                         NormVertexTabVector[i*10  ] = float(ValCol[j+1])*(fraction) + (1-fraction)*float(ValCol[(j-5)+1]);
 1751                         NormVertexTabVector[i*10+1] = float(ValCol[j+2])*(fraction) + (1-fraction)*float(ValCol[(j-5)+2]);;
 1752                         NormVertexTabVector[i*10+2] = float(ValCol[j+3])*(fraction) + (1-fraction)*float(ValCol[(j-5)+3]);
 1753                         NormVertexTabVector[i*10+3] = float(ValCol[(j)+4]);
 1754                         j = masterthread->VRgbtSize;
 1755                     }
 1756                 }
 1757                 else if(tmp == ValCol[j])
 1758                 {
 1759                     NormVertexTabVector[i*10  ] = float(ValCol[j+1]);
 1760                     NormVertexTabVector[i*10+1] = float(ValCol[j+2]);
 1761                     NormVertexTabVector[i*10+2] = float(ValCol[j+3]);
 1762                     NormVertexTabVector[i*10+3] = float(ValCol[j+4]);
 1763                     j = masterthread->VRgbtSize;
 1764                 }
 1765         }
 1766     }
 1767     else if(comp->ThereisRGBA[index] == true &&  comp->NoiseParam[comp->ParisoCurrentComponentIndex].NoiseType == 1)
 1768     {
 1769         uint idx=0;
 1770         for(uint i=0; i < comp->NbComponentsType.size()-1; i++)
 1771             idx+=comp->NbComponentsType[i];
 1772         for(uint i= comp->ParisoVertex[2*idx]; i < NbVertexTmp; i++)
 1773         {
 1774             if((i >= uint(comp->ParisoVertex[2*(cmpId+idx)])))
 1775             {
 1776                 K = cmpId;
 1777                 if((masterthread->componentsNumber-1)>cmpId)
 1778                 {
 1779                     cmpId++;
 1780                 }
 1781             }
 1782             val[0]= double(NormVertexTabVector[i*10 + 7]);
 1783             val[1]= double(NormVertexTabVector[i*10 + 8]);
 1784             val[2]= double(NormVertexTabVector[i*10 + 9]);
 1785             val[4]= double(K);
 1786             if(masterthread->gridnotnull)
 1787             {
 1788                 val[5] = double(i);
 1789                 val[6] = masterthread->x_Step[K];
 1790                 val[7] = masterthread->y_Step[K];
 1791                 val[8] = masterthread->z_Step[K];
 1792                 val[9] = double(masterthread->grid[K]);
 1793                 val[10] = masterthread->x_Sup[K];
 1794                 val[11] = masterthread->y_Sup[K];
 1795                 val[12] = masterthread->z_Sup[K];
 1796                 val[13] = masterthread->x_Inf[K];
 1797                 val[14] = masterthread->y_Inf[K];
 1798                 val[15] = masterthread->z_Inf[K];
 1799             }
 1800             else
 1801             {
 1802                 val[5] = double(i);
 1803                 val[6] = masterthread->x_Step[0];
 1804                 val[7] = masterthread->y_Step[0];
 1805                 val[8] = masterthread->z_Step[0];
 1806                 val[9] = double(masterthread->GridVal);
 1807                 val[10] = masterthread->x_Sup[0];
 1808                 val[11] = masterthread->y_Sup[0];
 1809                 val[12] = masterthread->z_Sup[0];
 1810                 val[13] = masterthread->x_Inf[0];
 1811                 val[14] = masterthread->y_Inf[0];
 1812                 val[15] = masterthread->z_Inf[0];
 1813             }
 1814             if(masterthread->Noise != "")
 1815                 tmp  = masterthread->NoiseParser->Eval(val);
 1816             else
 1817                 tmp =1.0;
 1818             val[0]= tmp*double(NormVertexTabVector[i*10+7]);
 1819             val[1]= tmp*double(NormVertexTabVector[i*10+8]);
 1820             val[2]= tmp*double(NormVertexTabVector[i*10+9]);
 1821             NormVertexTabVector[i*10  ] = float(masterthread->RgbtParser[0].Eval(val));
 1822             NormVertexTabVector[i*10+1] = float(masterthread->RgbtParser[1].Eval(val));
 1823             NormVertexTabVector[i*10+2] = float(masterthread->RgbtParser[2].Eval(val));
 1824             NormVertexTabVector[i*10+3] = float(masterthread->RgbtParser[3].Eval(val));
 1825         }
 1826     }
 1827     delete[] ValCol;
 1828 }
 1829 uint Iso3D::CNDCalculation(uint & NbTriangleIsoSurfaceTmp, struct ComponentInfos *comp)
 1830 {
 1831     uint idpx=0;
 1832     for(uint i=0; i < comp->NbComponentsType.size()-1; i++)
 1833         idpx+=comp->NbComponentsType[i];
 1834     uint startpoint=comp->ParisoVertex[2*idpx];
 1835     //In the case the parametric part of a Pariso object doesn't have a CND condition
 1836     int sz = (comp->ParisoCondition.size() ==
 1837               comp->NbComponentsType[comp->NbComponentsType.size()-1]) ? 0 : idpx;
 1838     if (masterthread->ParisoCondition == 1)
 1839     {
 1840         double vals[4];
 1841         std::vector<int> PointVerifyCond;
 1842         vals[3] = masterthread->stepMorph;
 1843         for(uint i= startpoint; i < NbVertexTmp; i++)
 1844         {
 1845             vals[0] = double(NormVertexTabVector[i*10+7]);
 1846             vals[1] = double(NormVertexTabVector[i*10+8]);
 1847             vals[2] = double(NormVertexTabVector[i*10+9]);
 1848             uint compid= CNDtoUse(i, comp);
 1849             if(comp->ParisoCondition[compid+sz])
 1850                 PointVerifyCond.push_back(8);
 1851             else
 1852                 PointVerifyCond.push_back(int(masterthread->ParisoConditionParser[compid].Eval(vals)));
 1853             if(PointVerifyCond[i-startpoint])
 1854             {
 1855                 NormVertexTabVector[i*10  ] = 0.1f;
 1856                 NormVertexTabVector[i*10+1] = 0.9f;
 1857                 NormVertexTabVector[i*10+2] = 0.0;
 1858                 NormVertexTabVector[i*10+3] = 1.0;
 1859             }
 1860             else
 1861             {
 1862                 NormVertexTabVector[i*10  ] = 0.9f;
 1863                 NormVertexTabVector[i*10+1] = 0.1f;
 1864                 NormVertexTabVector[i*10+2] = 0.9;
 1865                 NormVertexTabVector[i*10+3] = 1.0;
 1866             }
 1867         }
 1868         uint Aindex, Bindex, Cindex;
 1869         uint nbtriangle = NbTriangleIsoSurfaceTmp;
 1870         uint mdx=0;
 1871         for(uint id=0; id < comp->NbComponentsType.size()-1; id++)
 1872             mdx+=comp->NbComponentsType[id];
 1873         uint starttri = uint(comp->ParisoTriangle[2*mdx]/3);
 1874         std::vector<int> TypeIsoSurfaceTriangleListeCNDVector (NbTriangleIsoSurfaceTmp-starttri, 1);
 1875         for(uint i= starttri; i < nbtriangle; i++)
 1876         {
 1877             Aindex = IndexPolyTabVector[3*i    ];
 1878             Bindex = IndexPolyTabVector[3*i + 1];
 1879             Cindex = IndexPolyTabVector[3*i + 2];
 1880             int TypeTriangle = -1;
 1881             if((PointVerifyCond[Aindex-startpoint] == 8) ||
 1882                     (PointVerifyCond[Bindex-startpoint] == 8) ||
 1883                     (PointVerifyCond[Cindex-startpoint] == 8))
 1884             {
 1885                 TypeTriangle = 8;
 1886                 TypeIsoSurfaceTriangleListeCNDVector[i-starttri] = 8;
 1887             }
 1888             else if(PointVerifyCond[Aindex-startpoint] && !PointVerifyCond[Bindex-startpoint] && !PointVerifyCond[Cindex-startpoint])
 1889                 TypeTriangle = 0;
 1890             else if(!PointVerifyCond[Aindex-startpoint] && PointVerifyCond[Bindex-startpoint] && PointVerifyCond[Cindex-startpoint])
 1891                 TypeTriangle = 1;
 1892             else if(!PointVerifyCond[Aindex-startpoint] && PointVerifyCond[Bindex-startpoint] && !PointVerifyCond[Cindex-startpoint])
 1893                 TypeTriangle = 2;
 1894             else if(PointVerifyCond[Aindex-startpoint] && !PointVerifyCond[Bindex-startpoint] && PointVerifyCond[Cindex-startpoint])
 1895                 TypeTriangle = 3;
 1896             else if(!PointVerifyCond[Aindex-startpoint] && !PointVerifyCond[Bindex-startpoint] && PointVerifyCond[Cindex-startpoint])
 1897                 TypeTriangle = 4;
 1898             else if(PointVerifyCond[Aindex-startpoint] && PointVerifyCond[Bindex-startpoint] && !PointVerifyCond[Cindex-startpoint])
 1899                 TypeTriangle = 5;
 1900             else if(!PointVerifyCond[Aindex-startpoint] && !PointVerifyCond[Bindex-startpoint] && !PointVerifyCond[Cindex-startpoint])
 1901             {
 1902                 TypeTriangle = 6;
 1903                 TypeIsoSurfaceTriangleListeCNDVector[i-starttri] = -1;
 1904             }
 1905             else if(PointVerifyCond[Aindex-startpoint] && PointVerifyCond[Bindex-startpoint] && PointVerifyCond[Cindex-startpoint])
 1906             {
 1907                 TypeTriangle = 7;
 1908                 TypeIsoSurfaceTriangleListeCNDVector[i-starttri] = 1;
 1909             }
 1910 
 1911             if(TypeTriangle == 2 || TypeTriangle == 3)
 1912             {
 1913                 Aindex = IndexPolyTabVector[3*i+1];
 1914                 Bindex = IndexPolyTabVector[3*i+2];
 1915                 Cindex = IndexPolyTabVector[3*i  ];
 1916             }
 1917             else if(TypeTriangle == 4 || TypeTriangle == 5)
 1918             {
 1919                 Aindex = IndexPolyTabVector[3*i+2];
 1920                 Bindex = IndexPolyTabVector[3*i  ];
 1921                 Cindex = IndexPolyTabVector[3*i+1];
 1922             }
 1923 
 1924             double Bprime[4], Cprime[4], DiffX, DiffY, DiffZ;
 1925             int Alfa;
 1926             uint cnd = CNDtoUse(Aindex, comp);
 1927             if(TypeTriangle >=0 && TypeTriangle <= 5)
 1928             {
 1929                 /// Bprime
 1930                 Bprime[0] = double(NormVertexTabVector[10*Aindex+7]);
 1931                 Bprime[1] = double(NormVertexTabVector[10*Aindex+8]);
 1932                 Bprime[2] = double(NormVertexTabVector[10*Aindex+9]);
 1933                 Bprime[3] = masterthread->stepMorph;
 1934 
 1935                 DiffX = double(NormVertexTabVector[10*Bindex+7] - NormVertexTabVector[3+10*Aindex  + 4])/20.0;
 1936                 DiffY = double(NormVertexTabVector[10*Bindex+8] - NormVertexTabVector[3+10*Aindex+1+ 4])/20.0;
 1937                 DiffZ = double(NormVertexTabVector[10*Bindex+9] - NormVertexTabVector[3+10*Aindex+2+ 4])/20.0;
 1938                 Alfa = 0;
 1939                 if(TypeTriangle == 0 || TypeTriangle == 2 || TypeTriangle == 4)
 1940                 {
 1941                     while(masterthread->ParisoConditionParser[cnd].Eval(Bprime) == 1.0 && (Alfa < 20))
 1942                     {
 1943                         Bprime[0] += DiffX;
 1944                         Bprime[1] += DiffY;
 1945                         Bprime[2] += DiffZ;
 1946                         Alfa += 1;
 1947                     }
 1948                 }
 1949                 else
 1950                 {
 1951                     while(!(masterthread->ParisoConditionParser[cnd].Eval(Bprime) == 1.0) && (Alfa < 20))
 1952                     {
 1953                         Bprime[0] += DiffX;
 1954                         Bprime[1] += DiffY;
 1955                         Bprime[2] += DiffZ;
 1956                         Alfa += 1;
 1957                     }
 1958                 }
 1959 
 1960                 /// Cprime
 1961                 Cprime[0] = double(NormVertexTabVector[10*Aindex+7]);
 1962                 Cprime[1] = double(NormVertexTabVector[10*Aindex+8]);
 1963                 Cprime[2] = double(NormVertexTabVector[10*Aindex+9]);
 1964                 Cprime[3] = masterthread->stepMorph;
 1965 
 1966                 DiffX = double(NormVertexTabVector[3+10*Cindex  + 4] - NormVertexTabVector[3+10*Aindex  + 4])/20;
 1967                 DiffY = double(NormVertexTabVector[3+10*Cindex+1+ 4] - NormVertexTabVector[3+10*Aindex+1+ 4])/20;
 1968                 DiffZ = double(NormVertexTabVector[3+10*Cindex+2+ 4] - NormVertexTabVector[3+10*Aindex+2+ 4])/20;
 1969                 Alfa = 0;
 1970                 if(TypeTriangle == 0 || TypeTriangle == 2 || TypeTriangle == 4)
 1971                 {
 1972                     while(masterthread->ParisoConditionParser[cnd].Eval(Cprime) == 1.0 && (Alfa < 20))
 1973                     {
 1974                         Cprime[0] += DiffX;
 1975                         Cprime[1] += DiffY;
 1976                         Cprime[2] += DiffZ;
 1977                         Alfa += 1;
 1978                     }
 1979                 }
 1980                 else
 1981                 {
 1982                     while(!(masterthread->ParisoConditionParser[cnd].Eval(Cprime) == 1.0) && (Alfa < 20))
 1983                     {
 1984                         Cprime[0] += DiffX;
 1985                         Cprime[1] += DiffY;
 1986                         Cprime[2] += DiffZ;
 1987                         Alfa += 1;
 1988                     }
 1989                 }
 1990 
 1991                 //***********
 1992                 //Add points:
 1993                 //***********
 1994                 //Add Bprime:
 1995                 NormVertexTabVector.push_back(1.0);
 1996                 NormVertexTabVector.push_back(1.0);
 1997                 NormVertexTabVector.push_back(1.0);
 1998                 NormVertexTabVector.push_back(1.0);
 1999                 NormVertexTabVector.push_back(NormVertexTabVector[10*Bindex + 4]);
 2000                 NormVertexTabVector.push_back(NormVertexTabVector[10*Bindex + 5]);
 2001                 NormVertexTabVector.push_back(NormVertexTabVector[10*Bindex + 6]);
 2002                 NormVertexTabVector.push_back(float(Bprime[0]));
 2003                 NormVertexTabVector.push_back(float(Bprime[1]));
 2004                 NormVertexTabVector.push_back(float(Bprime[2]));
 2005 
 2006                 //Add Cprime:
 2007                 NormVertexTabVector.push_back(1.0);
 2008                 NormVertexTabVector.push_back(1.0);
 2009                 NormVertexTabVector.push_back(1.0);
 2010                 NormVertexTabVector.push_back(1.0);
 2011                 NormVertexTabVector.push_back(NormVertexTabVector[10*Cindex + 4]);
 2012                 NormVertexTabVector.push_back(NormVertexTabVector[10*Cindex + 5]);
 2013                 NormVertexTabVector.push_back(NormVertexTabVector[10*Cindex + 6]);
 2014                 NormVertexTabVector.push_back(float(Cprime[0]));
 2015                 NormVertexTabVector.push_back(float(Cprime[1]));
 2016                 NormVertexTabVector.push_back(float(Cprime[2]));
 2017 
 2018                 NbVertexTmp += 2;
 2019 
 2020                 //***********
 2021                 //Add triangles:
 2022                 //***********
 2023                 /// Add Three new triangles :
 2024                 uint IndexBprime = (NbVertexTmp-2);
 2025                 uint IndexCprime = (NbVertexTmp-1);
 2026 
 2027                 // The original triangle will be replaced by four other triangles:
 2028                 TypeIsoSurfaceTriangleListeCNDVector[i-starttri]=0;
 2029 
 2030                 /// (A, Bprime, Cprime)
 2031                 IndexPolyTabVector.push_back(Aindex);
 2032                 IndexPolyTabVector.push_back(IndexBprime);
 2033                 IndexPolyTabVector.push_back(IndexCprime);
 2034 
 2035                 (TypeTriangle == 0 || TypeTriangle == 2 || TypeTriangle == 4) ?
 2036                 TypeIsoSurfaceTriangleListeCNDVector.push_back(1) : TypeIsoSurfaceTriangleListeCNDVector.push_back(-1);
 2037                 NbTriangleIsoSurfaceTmp++;
 2038                 IndexPolyTabMinVector.push_back(3);
 2039                 IndexPolyTabMinVector.push_back(Aindex);
 2040                 IndexPolyTabMinVector.push_back(IndexBprime);
 2041                 IndexPolyTabMinVector.push_back(IndexCprime);
 2042 
 2043                 /// (Bprime, B, C)
 2044                 IndexPolyTabVector.push_back(IndexBprime);
 2045                 IndexPolyTabVector.push_back(Bindex);
 2046                 IndexPolyTabVector.push_back(Cindex);
 2047 
 2048                 (TypeTriangle == 0 || TypeTriangle == 2 || TypeTriangle == 4) ?
 2049                 TypeIsoSurfaceTriangleListeCNDVector.push_back(-1) : TypeIsoSurfaceTriangleListeCNDVector.push_back(1);
 2050                 NbTriangleIsoSurfaceTmp++;
 2051                 IndexPolyTabMinVector.push_back(3);
 2052                 IndexPolyTabMinVector.push_back(IndexBprime);
 2053                 IndexPolyTabMinVector.push_back(Bindex);
 2054                 IndexPolyTabMinVector.push_back(Cindex);
 2055 
 2056                 /// (Bprime, C, Cprime)
 2057                 IndexPolyTabVector.push_back(IndexBprime);
 2058                 IndexPolyTabVector.push_back(Cindex);
 2059                 IndexPolyTabVector.push_back(IndexCprime);
 2060 
 2061                 (TypeTriangle == 0 || TypeTriangle == 2 || TypeTriangle == 4) ?
 2062                 TypeIsoSurfaceTriangleListeCNDVector.push_back(-1) : TypeIsoSurfaceTriangleListeCNDVector.push_back(1);
 2063                 NbTriangleIsoSurfaceTmp++;
 2064                 IndexPolyTabMinVector.push_back(3);
 2065                 IndexPolyTabMinVector.push_back(IndexBprime);
 2066                 IndexPolyTabMinVector.push_back(Cindex);
 2067                 IndexPolyTabMinVector.push_back(IndexCprime);
 2068 
 2069                 /// (Bprime, Cprime) --> the border
 2070                 IndexPolyTabVector.push_back(IndexBprime);
 2071                 IndexPolyTabVector.push_back(IndexCprime);
 2072                 IndexPolyTabVector.push_back(IndexCprime);
 2073 
 2074                 TypeIsoSurfaceTriangleListeCNDVector.push_back(4); /// Type = 4-->Border
 2075                 NbTriangleIsoSurfaceTmp++;
 2076                 IndexPolyTabMinVector.push_back(3);
 2077                 IndexPolyTabMinVector.push_back(IndexBprime);
 2078                 IndexPolyTabMinVector.push_back(IndexCprime);
 2079                 IndexPolyTabMinVector.push_back(IndexCprime);
 2080             }
 2081         }
 2082 
 2083         //***********
 2084         //Reorganize the triangles index:
 2085         //***********
 2086         std::vector<uint>NewIndexPolyTabVector;
 2087         uint k, l, M, N;
 2088         k = l = M = N = 0;
 2089 
 2090         // In case we have a ParIso object, this will save the triangle arrangement for the parametric CND
 2091         for(uint i=0; i<starttri; i++)
 2092         {
 2093             NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i  ]);
 2094             NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+1]);
 2095             NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+2]);
 2096         }
 2097 
 2098         // Now we start sorting the triangles list. We have 3 cases
 2099 
 2100         for(uint i=starttri; i<NbTriangleIsoSurfaceTmp; i++)
 2101             if(TypeIsoSurfaceTriangleListeCNDVector[i-starttri] == 8)
 2102             {
 2103                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i  ]);
 2104                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+1]);
 2105                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+2]);
 2106                 N++;
 2107             }
 2108 
 2109         for(uint i=starttri; i<NbTriangleIsoSurfaceTmp; i++)
 2110             if(TypeIsoSurfaceTriangleListeCNDVector[i-starttri] == 1)
 2111             {
 2112                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i  ]);
 2113                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+1]);
 2114                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+2]);
 2115                 k++;
 2116             }
 2117 
 2118         for(uint i=starttri; i<NbTriangleIsoSurfaceTmp; i++)
 2119             if(TypeIsoSurfaceTriangleListeCNDVector[i-starttri] == -1)
 2120             {
 2121                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i  ]);
 2122                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+1]);
 2123                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+2]);
 2124                 l++;
 2125             }
 2126 
 2127         for(uint i=starttri; i<NbTriangleIsoSurfaceTmp; i++)
 2128             if(TypeIsoSurfaceTriangleListeCNDVector[i-starttri] == 4)
 2129             {
 2130                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i  ]);
 2131                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+1]);
 2132                 NewIndexPolyTabVector.push_back(IndexPolyTabVector[3*i+2]);
 2133                 M++;
 2134             }
 2135         //Copy the new index in the original one:
 2136         IndexPolyTabVector.clear();
 2137         NewIndexPolyTabVector.shrink_to_fit();
 2138         IndexPolyTabVector = NewIndexPolyTabVector;
 2139         NewIndexPolyTabVector.clear();
 2140         NewIndexPolyTabVector.shrink_to_fit();
 2141 
 2142         NbTriangleIsoSurfaceTmp = M + l + k + N;
 2143 
 2144         comp->NbTrianglesNoCND.push_back(N);
 2145         comp->NbTrianglesVerifyCND.push_back(k);
 2146         comp->NbTrianglesNotVerifyCND.push_back(l);
 2147         comp->NbTrianglesBorderCND.push_back(M);
 2148         comp->ThereisCND.push_back(true);
 2149 
 2150         PointVerifyCond.clear();
 2151         PointVerifyCond.shrink_to_fit();
 2152         TypeIsoSurfaceTriangleListeCNDVector.clear();
 2153         TypeIsoSurfaceTriangleListeCNDVector.shrink_to_fit();
 2154     }
 2155     else
 2156     {
 2157         comp->NbTrianglesNoCND.push_back(0);
 2158         comp->NbTrianglesVerifyCND.push_back(0);
 2159         comp->NbTrianglesNotVerifyCND.push_back(0);
 2160         comp->NbTrianglesBorderCND.push_back(0);
 2161         comp->ThereisCND.push_back(false);
 2162         for(uint i= startpoint; i < NbVertexTmp; i++)
 2163         {
 2164             NormVertexTabVector[i*10  ] = 0.5f;
 2165             NormVertexTabVector[i*10+1] = 0.6f;
 2166             NormVertexTabVector[i*10+2] = 0.8f;
 2167             NormVertexTabVector[i*10+3] = 1.0;
 2168         }
 2169     }
 2170 
 2171     return 1;
 2172 }
 2173 
 2174 Iso3D::~Iso3D()
 2175 {
 2176 }
 2177 
 2178 uint Iso3D::SetMiniMmeshStruct()
 2179 {
 2180     uint I, J, IJK, i, j, k, Index, lnew, iter, nbpl, iterpl;
 2181     uint maxgrscalemaxgr = masterthread->GridVal*masterthread->GridVal;
 2182 
 2183     lnew = 0;
 2184     NbPolyMin = 0;
 2185     /// Copy Index Polygons :
 2186     for(i=0; i+1 < masterthread->XYZgrid; i++)
 2187     {
 2188         I  = i*maxgrscalemaxgr;
 2189         for(j=0; j+1 < masterthread->XYZgrid; j++)
 2190         {
 2191             J  = I+j*masterthread->GridVal;
 2192             for(k=0; k+1 < masterthread->XYZgrid; k++)
 2193             {
 2194                 IJK = J+k;
 2195                 Index = GridVoxelVarPt[IJK].Signature;
 2196                 nbpl = triTable_min[Index][16];
 2197 
 2198                 if( nbpl == 1)
 2199                 {
 2200                     NbPolyMin += nbpl;
 2201                     IndexPolyTabMinVector.push_back(triTable_min[Index][17]);
 2202                     lnew++;
 2203                     for(iter = 0; iter < triTable_min[Index][17]; iter++)
 2204                         IndexPolyTabMinVector.push_back(GridVoxelVarPt[IJK].Edge_Points[triTable_min[Index][iter]]  + NbVertexTmp);
 2205                     lnew+=triTable_min[Index][17];
 2206                 }
 2207                 else if( nbpl == 2)
 2208                 {
 2209                     NbPolyMin += nbpl;
 2210                     /// First Poly:
 2211                     IndexPolyTabMinVector.push_back(triTable_min[Index][17]);
 2212                     lnew++;
 2213                     for(iter = 0; iter < triTable_min[Index][17]; iter++)
 2214                         IndexPolyTabMinVector.push_back(GridVoxelVarPt[IJK].Edge_Points[triTable_min[Index][iter]]  + NbVertexTmp);
 2215                     lnew+=triTable_min[Index][17];
 2216                     /// Second Poly:
 2217                     IndexPolyTabMinVector.push_back(triTable_min[Index][18]);
 2218                     lnew++;
 2219                     for(iter = triTable_min[Index][17]; iter < triTable_min[Index][17]+triTable_min[Index][18]; iter++)
 2220                         IndexPolyTabMinVector.push_back(GridVoxelVarPt[IJK].Edge_Points[triTable_min[Index][iter]]  + NbVertexTmp);
 2221                     lnew+=triTable_min[Index][17]+triTable_min[Index][18];
 2222                 }
 2223                 else if( nbpl > 2)
 2224                 {
 2225                     NbPolyMin += nbpl;
 2226                     /// In this case we have only Triangles (3 or 4):
 2227                     iter = 0;
 2228                     for(iterpl = 0; iterpl < nbpl; iterpl++)
 2229                     {
 2230                         IndexPolyTabMinVector.push_back(3);
 2231                         IndexPolyTabMinVector.push_back(GridVoxelVarPt[IJK].Edge_Points[triTable_min[Index][iter  ]]  + NbVertexTmp);
 2232                         IndexPolyTabMinVector.push_back(GridVoxelVarPt[IJK].Edge_Points[triTable_min[Index][iter+1]]  + NbVertexTmp);
 2233                         IndexPolyTabMinVector.push_back(GridVoxelVarPt[IJK].Edge_Points[triTable_min[Index][iter+2]]  + NbVertexTmp);
 2234                         lnew += 4;
 2235                         iter +=3;
 2236                     }
 2237                 }
 2238             } /// End of for(k=0;
 2239         } /// End of for(j=0;
 2240     } /// End of for(i=0;
 2241     return 1;
 2242 }
 2243 
 2244 uint Iso3D::ConstructIsoSurface()
 2245 {
 2246     uint Index, IndexFirstPoint, IndexSeconPoint, IndexThirdPoint;
 2247     uint I, J, IJK;
 2248     uint maxgrscalemaxgr = masterthread->GridVal*masterthread->GridVal;
 2249 
 2250     NbTriangleIsoSurface = 0;
 2251     for(uint i=0; i+1 < masterthread->XYZgrid; i++)
 2252     {
 2253         I   = i*maxgrscalemaxgr;
 2254         for(uint j=0; j+1 < masterthread->XYZgrid; j++)
 2255         {
 2256             J   = I+j*masterthread->GridVal;
 2257             for(uint k=0; k+1 < masterthread->XYZgrid; k++)
 2258             {
 2259                 IJK = J+k;
 2260                 Index = GridVoxelVarPt[IJK].Signature;
 2261                 for (uint l = 0; triTable[Index][l] != 111; l += 3)
 2262                 {
 2263                     IndexFirstPoint = GridVoxelVarPt[IJK].Edge_Points[triTable[Index][l  ]];
 2264                     IndexSeconPoint = GridVoxelVarPt[IJK].Edge_Points[triTable[Index][l+1]];
 2265                     IndexThirdPoint = GridVoxelVarPt[IJK].Edge_Points[triTable[Index][l+2]];
 2266                     IndexPolyTabVector.push_back(IndexFirstPoint + NbVertexTmp);
 2267                     IndexPolyTabVector.push_back(IndexSeconPoint + NbVertexTmp);
 2268                     IndexPolyTabVector.push_back(IndexThirdPoint + NbVertexTmp);
 2269                     NbTriangleIsoSurface++;
 2270                 }
 2271             }
 2272         }
 2273     }
 2274     return 1;
 2275 }
 2276 
 2277 uint Iso3D::PointEdgeComputation(uint isoindex)
 2278 {
 2279     uint i_Start, i_End, j_Start, j_End, k_Start, k_End, i, j, k;
 2280     double vals[7], IsoValue_1, IsoValue_2, rapport;
 2281     double factor;
 2282     uint maxgridval=masterthread->GridVal;
 2283     uint maxgrscalemaxgr = maxgridval*maxgridval;
 2284     uint I, J, JK, IJK, IPLUSONEJK, IMINUSONEJK,
 2285          IJPLUSONEK, IJMINUSONEK, IMINUSONEJMINUSONEK, IsoValue=0, NbPointIsoMap_local=0;
 2286     /// We have to compute the edges for the Grid limits ie: i=0, j=0 and k=0
 2287 
 2288     i_Start = 1/*+masterthread->iStart*/;
 2289     j_Start = 1;
 2290     k_Start = 1;
 2291 
 2292     i_End = masterthread->XYZgrid-1;
 2293     j_End = masterthread->XYZgrid-1;
 2294     k_End = masterthread->XYZgrid-1;
 2295     /// The code is doubled to eliminate conditions tests
 2296 
 2297     for(i = i_Start; i < i_End; i++)
 2298     {
 2299         I = i*maxgrscalemaxgr;
 2300         for(j = j_Start; j < j_End; j++)
 2301         {
 2302             J = I + j*maxgridval;
 2303             for(k = k_Start; k < k_End; k++)
 2304             {
 2305                 IJK                 = J+k;
 2306                 IPLUSONEJK          = IJK + maxgrscalemaxgr;
 2307                 IMINUSONEJK         = IJK - maxgrscalemaxgr;
 2308                 IJPLUSONEK          = IJK + maxgridval;
 2309                 IJMINUSONEK         = IJK - maxgridval;
 2310                 IMINUSONEJMINUSONEK = IMINUSONEJK - maxgridval;
 2311 
 2312                 IsoValue_1 = Results[IJK];
 2313                 /// First Case P(i+1)(j)(k)
 2314 
 2315                 IsoValue_2 = Results[IPLUSONEJK];
 2316                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0 )
 2317                 {
 2318                     // Edge Point computation and  save in IsoPointMap
 2319                     factor = (IsoValue - IsoValue_1)/rapport;
 2320 
 2321                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i] - factor * masterthread->x_Step[isoindex];
 2322                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2323                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2324                     ///===========================================================///
 2325                     for(int iter=0; iter<7; iter++)
 2326                         NormVertexTabVector.push_back(1.0);
 2327                     NormVertexTabVector.push_back(float(vals[0]));
 2328                     NormVertexTabVector.push_back(float(vals[1]));
 2329                     NormVertexTabVector.push_back(float(vals[2]));
 2330 
 2331                     // save The reference to this point
 2332                     GridVoxelVarPt[IJK].Edge_Points[0] = NbPointIsoMap_local;
 2333                     GridVoxelVarPt[IJK].NbEdgePoint += 1;
 2334 
 2335                     // The same Point is used in three other Voxels
 2336                     GridVoxelVarPt[IJMINUSONEK].Edge_Points[4] = NbPointIsoMap_local;
 2337                     GridVoxelVarPt[IJMINUSONEK].NbEdgePoint += 1;
 2338 
 2339                     GridVoxelVarPt[IJK-1].Edge_Points[2]   = NbPointIsoMap_local;
 2340                     GridVoxelVarPt[IJK-1].NbEdgePoint += 1;
 2341 
 2342                     GridVoxelVarPt[IJMINUSONEK-1].Edge_Points[6] = NbPointIsoMap_local;
 2343                     GridVoxelVarPt[IJMINUSONEK-1].NbEdgePoint += 1;
 2344 
 2345                     NbPointIsoMap_local++;
 2346                 }
 2347                 ///+++++++++++++++++++++++++++++++++++++++++
 2348                 /// Second Case P(i)(j+1)(k)
 2349 
 2350                 IsoValue_2 = Results[IJPLUSONEK];
 2351                 // Edge Point computation and  save in IsoPointMap
 2352                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2353                 {
 2354                     factor = (IsoValue - IsoValue_1)/rapport;
 2355 
 2356                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2357                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j] - factor * masterthread->y_Step[isoindex];
 2358                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2359                     for(int iter=0; iter<7; iter++)
 2360                         NormVertexTabVector.push_back(1.0);
 2361                     NormVertexTabVector.push_back(float(vals[0]));
 2362                     NormVertexTabVector.push_back(float(vals[1]));
 2363                     NormVertexTabVector.push_back(float(vals[2]));
 2364 
 2365                     // save The reference to this point
 2366                     GridVoxelVarPt[IJK].Edge_Points[8] = NbPointIsoMap_local;
 2367                     GridVoxelVarPt[IJK].NbEdgePoint += 1;
 2368 
 2369                     // The same Point is used in three other Voxels
 2370                     GridVoxelVarPt[IMINUSONEJK].Edge_Points[9]    = NbPointIsoMap_local;
 2371                     GridVoxelVarPt[IMINUSONEJK].NbEdgePoint += 1;
 2372 
 2373                     GridVoxelVarPt[IJK-1].Edge_Points[11]   = NbPointIsoMap_local;
 2374                     GridVoxelVarPt[IJK-1].NbEdgePoint += 1;
 2375 
 2376                     GridVoxelVarPt[IMINUSONEJK-1].Edge_Points[10] = NbPointIsoMap_local;
 2377                     GridVoxelVarPt[IMINUSONEJK-1].NbEdgePoint += 1;
 2378 
 2379                     NbPointIsoMap_local++;
 2380                 }
 2381 
 2382                 // Third Case P(i)(j)(k+1)
 2383 
 2384                 IsoValue_2 = Results[IJK+1];
 2385                 // Edge Point computation and  save in IsoPointMap
 2386                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2387                 {
 2388                     factor = (IsoValue - IsoValue_1)/rapport;
 2389 
 2390                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2391                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2392                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k] - factor * masterthread->z_Step[isoindex];
 2393                     for(int iter=0; iter<7; iter++)
 2394                         NormVertexTabVector.push_back(1.0);
 2395                     NormVertexTabVector.push_back(float(vals[0]));
 2396                     NormVertexTabVector.push_back(float(vals[1]));
 2397                     NormVertexTabVector.push_back(float(vals[2]));
 2398                     // save The reference to this point
 2399                     GridVoxelVarPt[IJK].Edge_Points[3] = NbPointIsoMap_local;
 2400                     GridVoxelVarPt[IJK].NbEdgePoint += 1;
 2401                     // The same Point is used in three other Voxels
 2402                     GridVoxelVarPt[IMINUSONEJK].Edge_Points[1]   = NbPointIsoMap_local;
 2403                     GridVoxelVarPt[IMINUSONEJK].NbEdgePoint += 1;
 2404                     GridVoxelVarPt[IJMINUSONEK].Edge_Points[7]   = NbPointIsoMap_local;
 2405                     GridVoxelVarPt[IJMINUSONEK].NbEdgePoint += 1;
 2406                     GridVoxelVarPt[IMINUSONEJMINUSONEK].Edge_Points[5] = NbPointIsoMap_local;
 2407                     GridVoxelVarPt[IMINUSONEJMINUSONEK].NbEdgePoint += 1;
 2408                     NbPointIsoMap_local++;
 2409                 }
 2410             }
 2411         }
 2412     }
 2413 
 2414     /// Now we have to compute the Grid's limits...
 2415     /// The code is quite big but this is much more easy to compute
 2416 
 2417     /// 1) First case : i =0;
 2418 
 2419     i =0;
 2420     for(j=0; j < masterthread->XYZgrid; j++)
 2421     {
 2422         J = j*maxgridval;
 2423         for(k=0; k < masterthread->XYZgrid; k++)
 2424         {
 2425             JK = J +k;
 2426 
 2427             IsoValue_1 = Results[JK];
 2428 
 2429             IsoValue_2 = Results[maxgrscalemaxgr+JK];
 2430             if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2431             {
 2432                 // Edge Point computation and  save in IsoPointMap
 2433                 factor = (IsoValue - IsoValue_1)/rapport;
 2434 
 2435                 vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i] - factor * masterthread->x_Step[isoindex];
 2436                 vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2437                 vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2438                 for(int iter=0; iter<7; iter++)
 2439                     NormVertexTabVector.push_back(1.0);
 2440                 NormVertexTabVector.push_back(float(vals[0]));
 2441                 NormVertexTabVector.push_back(float(vals[1]));
 2442                 NormVertexTabVector.push_back(float(vals[2]));
 2443 
 2444                 // save The reference to this point
 2445                 GridVoxelVarPt[JK].Edge_Points[0] = NbPointIsoMap_local;
 2446                 GridVoxelVarPt[JK].NbEdgePoint += 1;
 2447 
 2448                 // The same Point is used in three other Voxels
 2449                 if( j != 0)
 2450                 {
 2451                     GridVoxelVarPt[JK-maxgridval].Edge_Points[4] = NbPointIsoMap_local;
 2452                     GridVoxelVarPt[JK-maxgridval].NbEdgePoint += 1;
 2453                 }
 2454                 if ( k != 0 )
 2455                 {
 2456                     GridVoxelVarPt[JK-1].Edge_Points[2]   = NbPointIsoMap_local;
 2457                     GridVoxelVarPt[JK-1].NbEdgePoint += 1;
 2458                 }
 2459                 if( j != 0 && k != 0)
 2460                 {
 2461                     GridVoxelVarPt[JK-maxgridval-1].Edge_Points[6] = NbPointIsoMap_local;
 2462                     GridVoxelVarPt[JK-maxgridval-1].NbEdgePoint += 1;
 2463                 }
 2464 
 2465                 NbPointIsoMap_local++;
 2466             }
 2467 
 2468             // Second Case P(0)(j+1)(k)
 2469             if ( j != (masterthread->XYZgrid -1))
 2470             {
 2471                 IsoValue_2 = Results[JK+maxgridval];
 2472                 // Edge Point computation and  save in IsoPointMap
 2473                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2474                 {
 2475                     factor = (IsoValue - IsoValue_1)/rapport;
 2476 
 2477                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2478                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j] - factor *masterthread->y_Step[isoindex];
 2479                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2480                     for(int iter=0; iter<7; iter++)
 2481                         NormVertexTabVector.push_back(1.0);
 2482                     NormVertexTabVector.push_back(float(vals[0]));
 2483                     NormVertexTabVector.push_back(float(vals[1]));
 2484                     NormVertexTabVector.push_back(float(vals[2]));
 2485                     // save The reference to this point
 2486                     GridVoxelVarPt[JK].Edge_Points[8] = NbPointIsoMap_local;
 2487                     GridVoxelVarPt[JK].NbEdgePoint += 1;
 2488 
 2489                     // The same Point is used in three other Voxels
 2490                     if ( k != 0)
 2491                     {
 2492                         GridVoxelVarPt[JK-1].Edge_Points[11]   = NbPointIsoMap_local;
 2493                         GridVoxelVarPt[JK-1].NbEdgePoint += 1;
 2494                     }
 2495 
 2496                     NbPointIsoMap_local++;
 2497                 }
 2498             } /// If ( j != nb_colon -1) ...
 2499 
 2500             // Third Case P(0)(j)(k+1)
 2501             if ( k != (masterthread->XYZgrid-1))
 2502             {
 2503                 IsoValue_2 = Results[JK+1];
 2504                 // Edge Point computation and  save in IsoPointMap
 2505                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2506                 {
 2507                     factor = (IsoValue - IsoValue_1)/rapport;
 2508 
 2509                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2510                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2511                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k] - factor * masterthread->z_Step[isoindex];
 2512                     for(int iter=0; iter<7; iter++)
 2513                         NormVertexTabVector.push_back(1.0);
 2514                     NormVertexTabVector.push_back(float(vals[0]));
 2515                     NormVertexTabVector.push_back(float(vals[1]));
 2516                     NormVertexTabVector.push_back(float(vals[2]));
 2517 
 2518                     // save The reference to this point
 2519                     GridVoxelVarPt[JK].Edge_Points[3] = NbPointIsoMap_local;
 2520                     GridVoxelVarPt[JK].NbEdgePoint += 1;
 2521 
 2522                     // The same Point is used in three other Voxels
 2523                     if ( j != 0)
 2524                     {
 2525                         GridVoxelVarPt[JK-maxgridval].Edge_Points[7]   =NbPointIsoMap_local;
 2526                         GridVoxelVarPt[JK-maxgridval].NbEdgePoint += 1;
 2527                     }
 2528 
 2529                     NbPointIsoMap_local++;
 2530 
 2531                 }
 2532             } /// End of ( if ( k != nb_depth -1)....
 2533         }
 2534     }
 2535 
 2536 
 2537 
 2538     /// 2) Case i = nb_ligne-1
 2539 
 2540     i = masterthread->XYZgrid-1;
 2541     I = i*maxgrscalemaxgr;
 2542     for(j=0; j < masterthread->XYZgrid; j++)
 2543     {
 2544         J = I + j*maxgridval;
 2545         for(k=0; k < masterthread->XYZgrid; k++)
 2546         {
 2547             JK = J + k;
 2548 
 2549             IsoValue_1 = Results[JK];
 2550 
 2551 
 2552             // Second Case P(i)(j+1)(k)
 2553             if ( j != (masterthread->XYZgrid -1))
 2554             {
 2555                 IsoValue_2 = Results[JK+maxgridval];
 2556                 // Edge Point computation and  save in IsoPointMap
 2557                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2558                 {
 2559                     factor = (IsoValue - IsoValue_1)/rapport;
 2560 
 2561                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2562                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j] - factor *masterthread->y_Step[isoindex];
 2563                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2564                     for(int iter=0; iter<7; iter++)
 2565                         NormVertexTabVector.push_back(1.0);
 2566                     NormVertexTabVector.push_back(float(vals[0]));
 2567                     NormVertexTabVector.push_back(float(vals[1]));
 2568                     NormVertexTabVector.push_back(float(vals[2]));
 2569 
 2570                     // save The reference to this point
 2571                     GridVoxelVarPt[JK].Edge_Points[8] = NbPointIsoMap_local;
 2572                     GridVoxelVarPt[JK].NbEdgePoint += 1;
 2573 
 2574                     // The same Point is used in three other Voxels
 2575                     if(i != 0)
 2576                     {
 2577                         GridVoxelVarPt[JK-maxgrscalemaxgr].Edge_Points[9]    =NbPointIsoMap_local;
 2578                         GridVoxelVarPt[JK-maxgrscalemaxgr].NbEdgePoint += 1;
 2579                     }
 2580                     if(k != 0)
 2581                     {
 2582                         GridVoxelVarPt[JK-1].Edge_Points[11]   = NbPointIsoMap_local;
 2583                         GridVoxelVarPt[JK-1].NbEdgePoint += 1;
 2584                     }
 2585                     if(i != 0 && k != 0)
 2586                     {
 2587                         GridVoxelVarPt[JK-maxgrscalemaxgr-1].Edge_Points[10] = NbPointIsoMap_local;
 2588                         GridVoxelVarPt[JK-maxgrscalemaxgr-1].NbEdgePoint += 1;
 2589                     }
 2590 
 2591                     NbPointIsoMap_local++;
 2592 
 2593                 }
 2594             } /// End of if (j != nb_colon -1)...
 2595 
 2596             // Third Case P(i)(j)(k+1)
 2597             if ( k != (masterthread->XYZgrid -1))
 2598             {
 2599                 IsoValue_2 = Results[JK+1];
 2600                 // Edge Point computation and  save in IsoPointMap
 2601                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2602                 {
 2603                     factor = (IsoValue - IsoValue_1)/rapport;
 2604 
 2605                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2606                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2607                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k] - factor * masterthread->z_Step[isoindex];
 2608                     for(int iter=0; iter<7; iter++)
 2609                         NormVertexTabVector.push_back(1.0);
 2610                     NormVertexTabVector.push_back(float(vals[0]));
 2611                     NormVertexTabVector.push_back(float(vals[1]));
 2612                     NormVertexTabVector.push_back(float(vals[2]));
 2613 
 2614                     // save The reference to this point
 2615                     GridVoxelVarPt[JK].Edge_Points[3] = NbPointIsoMap_local;
 2616                     GridVoxelVarPt[JK].NbEdgePoint += 1;
 2617 
 2618                     // The same Point is used in three other Voxels
 2619                     if( i != 0)
 2620                     {
 2621                         GridVoxelVarPt[JK-maxgrscalemaxgr].Edge_Points[1]   = NbPointIsoMap_local;
 2622                         GridVoxelVarPt[JK-maxgrscalemaxgr].NbEdgePoint += 1;
 2623                     }
 2624                     if(j != 0)
 2625                     {
 2626                         GridVoxelVarPt[JK-maxgridval].Edge_Points[7]   = NbPointIsoMap_local;
 2627                         GridVoxelVarPt[JK-maxgridval].NbEdgePoint += 1;
 2628                     }
 2629                     if( i !=0 && j != 0)
 2630                     {
 2631                         GridVoxelVarPt[JK-maxgrscalemaxgr-maxgridval].Edge_Points[5] = NbPointIsoMap_local;
 2632                         GridVoxelVarPt[JK-maxgrscalemaxgr-maxgridval].NbEdgePoint += 1;
 2633                     }
 2634 
 2635                     NbPointIsoMap_local++;
 2636 
 2637                 }
 2638             } /// End of if ( k != nb_depth -1)...
 2639         }
 2640     }
 2641 
 2642     /// 3) Case j = 0
 2643     j= 0;
 2644 
 2645     for(i=0; i < masterthread->XYZgrid; i++)
 2646         for(k=0; k < masterthread->XYZgrid; k++)
 2647         {
 2648 
 2649             IsoValue_1 = Results[i*maxgrscalemaxgr+k];
 2650 
 2651             // First Case P(i+1)(j)(k)
 2652             if( i != (masterthread->XYZgrid -1))
 2653             {
 2654                 IsoValue_2 = Results[(i+1)*maxgrscalemaxgr+k];
 2655                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2656                 {
 2657 
 2658                     // Edge Point computation and  save in IsoPointMap
 2659                     factor = (IsoValue - IsoValue_1)/rapport;
 2660 
 2661                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i] - factor * masterthread->x_Step[isoindex];
 2662                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2663                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2664                     for(int iter=0; iter<7; iter++)
 2665                         NormVertexTabVector.push_back(1.0);
 2666                     NormVertexTabVector.push_back(float(vals[0]));
 2667                     NormVertexTabVector.push_back(float(vals[1]));
 2668                     NormVertexTabVector.push_back(float(vals[2]));
 2669 
 2670                     // save The reference to this point
 2671                     GridVoxelVarPt[i*maxgrscalemaxgr+k].Edge_Points[0] = NbPointIsoMap_local;
 2672                     GridVoxelVarPt[i*maxgrscalemaxgr+k].NbEdgePoint += 1;
 2673 
 2674                     // The same Point is used in three other Voxels
 2675                     if( k!=0)
 2676                     {
 2677                         GridVoxelVarPt[i*maxgrscalemaxgr+k-1].Edge_Points[2]   = NbPointIsoMap_local;
 2678                         GridVoxelVarPt[i*maxgrscalemaxgr+k-1].NbEdgePoint += 1;
 2679                     }
 2680 
 2681                     NbPointIsoMap_local++;
 2682                 }
 2683             } /// End of if ( i != nb_ligne -1)...
 2684 
 2685             // Second Case P(i)(j+1)(k)
 2686 
 2687             IsoValue_2 = Results[i*maxgrscalemaxgr+maxgridval+k];
 2688             // Edge Point computation and  save in IsoPointMap
 2689             if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2690             {
 2691                 factor = (IsoValue - IsoValue_1)/rapport;
 2692 
 2693                 vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2694                 vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j] - factor * masterthread->y_Step[isoindex];
 2695                 vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2696                 for(int iter=0; iter<7; iter++)
 2697                     NormVertexTabVector.push_back(1.0);
 2698                 NormVertexTabVector.push_back(float(vals[0]));
 2699                 NormVertexTabVector.push_back(float(vals[1]));
 2700                 NormVertexTabVector.push_back(float(vals[2]));
 2701 
 2702                 // save The reference to this point
 2703                 GridVoxelVarPt[i*maxgrscalemaxgr+k].Edge_Points[8] = NbPointIsoMap_local;
 2704                 GridVoxelVarPt[i*maxgrscalemaxgr+k].NbEdgePoint += 1;
 2705 
 2706                 // The same Point is used in three other Voxels
 2707                 if( i !=0 )
 2708                 {
 2709                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+k].Edge_Points[9]    = NbPointIsoMap_local;
 2710                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+k].NbEdgePoint += 1;
 2711                 }
 2712                 if(k !=0)
 2713                 {
 2714                     GridVoxelVarPt[i*maxgrscalemaxgr+k-1].Edge_Points[11]   = NbPointIsoMap_local;
 2715                     GridVoxelVarPt[i*maxgrscalemaxgr+k-1].NbEdgePoint += 1;
 2716                 }
 2717                 if( i !=0 && k !=0)
 2718                 {
 2719                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+k-1].Edge_Points[10] = NbPointIsoMap_local;
 2720                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+k-1].NbEdgePoint += 1;
 2721                 }
 2722                 NbPointIsoMap_local++;
 2723             }
 2724 
 2725 
 2726             // Third Case P(i)(j)(k+1)
 2727             if(k != (masterthread->XYZgrid -1))
 2728             {
 2729                 IsoValue_2 = Results[i*maxgrscalemaxgr+k+1];
 2730                 // Edge Point computation and  save in IsoPointMap
 2731                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2732                 {
 2733                     factor = (IsoValue - IsoValue_1)/rapport;
 2734 
 2735                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2736                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2737                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k] - factor * masterthread->z_Step[isoindex];
 2738                     for(int iter=0; iter<7; iter++)
 2739                         NormVertexTabVector.push_back(1.0);
 2740                     NormVertexTabVector.push_back(float(vals[0]));
 2741                     NormVertexTabVector.push_back(float(vals[1]));
 2742                     NormVertexTabVector.push_back(float(vals[2]));
 2743 
 2744                     // save The reference to this point
 2745                     GridVoxelVarPt[i*maxgrscalemaxgr+k].Edge_Points[3] = NbPointIsoMap_local;
 2746                     GridVoxelVarPt[i*maxgrscalemaxgr+k].NbEdgePoint += 1;
 2747 
 2748                     // The same Point is used in three other Voxels
 2749                     if( i != 0)
 2750                     {
 2751                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+k].Edge_Points[1]   = NbPointIsoMap_local;
 2752                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+k].NbEdgePoint += 1;
 2753                     }
 2754 
 2755                     NbPointIsoMap_local++;
 2756 
 2757                 }
 2758             } /// End of if(k != (nb_depth -1))...
 2759         }
 2760     /// 4) Case j = nb_colon -1
 2761     j = masterthread->XYZgrid-1;
 2762 
 2763     for(i=0; i < masterthread->XYZgrid; i++)
 2764         for(k=0; k < masterthread->XYZgrid; k++)
 2765         {
 2766             IsoValue_1 = Results[i*maxgrscalemaxgr+j*maxgridval+k];
 2767             // First Case P(i+1)(j)(k)
 2768             if( i != (masterthread->XYZgrid-1))
 2769             {
 2770                 IsoValue_2 = Results[(i+1)*maxgrscalemaxgr+j*maxgridval+k];
 2771                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2772                 {
 2773 
 2774                     // Edge Point computation and  save in IsoPointMap
 2775                     factor = (IsoValue - IsoValue_1)/rapport;
 2776 
 2777                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i] - factor * masterthread->x_Step[isoindex];
 2778                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2779                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2780                     for(int iter=0; iter<7; iter++)
 2781                         NormVertexTabVector.push_back(1.0);
 2782                     NormVertexTabVector.push_back(float(vals[0]));
 2783                     NormVertexTabVector.push_back(float(vals[1]));
 2784                     NormVertexTabVector.push_back(float(vals[2]));
 2785                     // save The reference to this point
 2786                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].Edge_Points[0] = NbPointIsoMap_local;
 2787                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].NbEdgePoint += 1;
 2788                     // The same Point is used in three other Voxels
 2789                     if( j != 0)
 2790                     {
 2791                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k].Edge_Points[4] = NbPointIsoMap_local;
 2792                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k].NbEdgePoint += 1;
 2793                     }
 2794                     if(k != 0)
 2795                     {
 2796                         GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k-1].Edge_Points[2]   = NbPointIsoMap_local;
 2797                         GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k-1].NbEdgePoint += 1;
 2798                     }
 2799                     if(j != 0 && k != 0)
 2800                     {
 2801                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k-1].Edge_Points[6] = NbPointIsoMap_local;
 2802                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k-1].NbEdgePoint += 1;
 2803                     }
 2804                     NbPointIsoMap_local++;
 2805                 }
 2806             } /// End of if( i != (nb_ligne-1))...
 2807 
 2808             // Third Case P(i)(j)(k+1)
 2809             if( k != (masterthread->XYZgrid -1))
 2810             {
 2811                 IsoValue_2 = Results[i*maxgrscalemaxgr+j*maxgridval+k+1];
 2812                 // Edge Point computation and  save in IsoPointMap
 2813                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2814                 {
 2815                     factor = (IsoValue - IsoValue_1)/rapport;
 2816 
 2817                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2818                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2819                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k] - factor * masterthread->z_Step[isoindex];
 2820                     for(int iter=0; iter<7; iter++)
 2821                         NormVertexTabVector.push_back(1.0);
 2822                     NormVertexTabVector.push_back(float(vals[0]));
 2823                     NormVertexTabVector.push_back(float(vals[1]));
 2824                     NormVertexTabVector.push_back(float(vals[2]));
 2825                     // save The reference to this point
 2826                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].Edge_Points[3] = NbPointIsoMap_local;
 2827                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].NbEdgePoint += 1;
 2828                     // The same Point is used in three other Voxels
 2829                     if(i !=0)
 2830                     {
 2831                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval+k].Edge_Points[1]   = NbPointIsoMap_local;
 2832                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval+k].NbEdgePoint += 1;
 2833                     }
 2834                     if( j!=0)
 2835                     {
 2836                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k].Edge_Points[7]   = NbPointIsoMap_local;
 2837                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k].NbEdgePoint += 1;
 2838                     }
 2839                     if(i != 0 && j !=0)
 2840                     {
 2841                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+(j-1)*maxgridval+k].Edge_Points[5] = NbPointIsoMap_local;
 2842                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+(j-1)*maxgridval+k].NbEdgePoint += 1;
 2843                     }
 2844                     NbPointIsoMap_local++;
 2845                 }
 2846             } /// End of if (k != nb_depth)...
 2847         }
 2848 
 2849     /// 5) Case k = 0
 2850 
 2851     k =0;
 2852 
 2853     for(i=0; i < masterthread->XYZgrid; i++)
 2854         for(j=0; j < masterthread->XYZgrid; j++)
 2855         {
 2856             IsoValue_1 = Results[i*maxgrscalemaxgr+j*maxgridval];
 2857             // First Case P(i+1)(j)(k)
 2858             if(i != (masterthread->XYZgrid -1))
 2859             {
 2860                 IsoValue_2 = Results[(i+1)*maxgrscalemaxgr+j*maxgridval];
 2861                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2862                 {
 2863                     // Edge Point computation and  save in IsoPointMap
 2864                     factor = (IsoValue - IsoValue_1)/rapport;
 2865 
 2866                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i] - factor * masterthread->x_Step[isoindex];
 2867                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2868                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2869                     for(int iter=0; iter<7; iter++)
 2870                         NormVertexTabVector.push_back(1.0);
 2871                     NormVertexTabVector.push_back(float(vals[0]));
 2872                     NormVertexTabVector.push_back(float(vals[1]));
 2873                     NormVertexTabVector.push_back(float(vals[2]));
 2874                     // save The reference to this point
 2875                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval].Edge_Points[0] = NbPointIsoMap_local;
 2876                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval].NbEdgePoint += 1;
 2877                     // The same Point is used in one other Voxel
 2878                     if(j !=0)
 2879                     {
 2880                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval].Edge_Points[4] = NbPointIsoMap_local;
 2881                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval].NbEdgePoint += 1;
 2882                     }
 2883                     NbPointIsoMap_local++;
 2884                 }
 2885             } /// End of if(i != (nb_ligne -1))
 2886 
 2887             // Second Case P(i)(j+1)(k)
 2888             if(j != masterthread->XYZgrid -1)
 2889             {
 2890                 IsoValue_2 = Results[i*maxgrscalemaxgr+(j+1)*maxgridval];
 2891                 // Edge Point computation and  save in IsoPointMap
 2892                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2893                 {
 2894                     factor = (IsoValue - IsoValue_1)/rapport;
 2895 
 2896                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2897                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j] - factor * masterthread->y_Step[isoindex];
 2898                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2899                     for(int iter=0; iter<7; iter++)
 2900                         NormVertexTabVector.push_back(1.0);
 2901                     NormVertexTabVector.push_back(float(vals[0]));
 2902                     NormVertexTabVector.push_back(float(vals[1]));
 2903                     NormVertexTabVector.push_back(float(vals[2]));
 2904                     // save The reference to this point
 2905                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval].Edge_Points[8] = NbPointIsoMap_local;
 2906                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval].NbEdgePoint += 1;
 2907                     // The same Point is used in one other Voxels
 2908                     if ( i !=0 )
 2909                     {
 2910                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval].Edge_Points[9]    = NbPointIsoMap_local;
 2911                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval].NbEdgePoint += 1;
 2912                     }
 2913                     NbPointIsoMap_local++;
 2914                 }
 2915             }/// End of if(j != nb_colon -1)...
 2916             // Third Case P(i)(j)(k+1)
 2917             IsoValue_2 = Results[i*maxgrscalemaxgr+j*maxgridval+1];
 2918             // Edge Point computation and  save in IsoPointMap
 2919             if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2920             {
 2921                 factor = (IsoValue - IsoValue_1)/rapport;
 2922 
 2923                 vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 2924                 vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2925                 vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k] - factor * masterthread->z_Step[isoindex];
 2926                 for(int iter=0; iter<7; iter++)
 2927                     NormVertexTabVector.push_back(1.0);
 2928                 NormVertexTabVector.push_back(float(vals[0]));
 2929                 NormVertexTabVector.push_back(float(vals[1]));
 2930                 NormVertexTabVector.push_back(float(vals[2]));
 2931                 // save The reference to this point
 2932                 GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval].Edge_Points[3] = NbPointIsoMap_local;
 2933                 GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval].NbEdgePoint += 1;
 2934 
 2935                 // The same Point is used in three other Voxels
 2936                 if ( i !=0 )
 2937                 {
 2938                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval].Edge_Points[1]   = NbPointIsoMap_local;
 2939                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval].NbEdgePoint += 1;
 2940                 }
 2941                 if (j != 0 )
 2942                 {
 2943                     GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval].Edge_Points[7]   = NbPointIsoMap_local;
 2944                     GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval].NbEdgePoint += 1;
 2945                 }
 2946                 if ( i !=0 && j != 0 )
 2947                 {
 2948                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+(j-1)*maxgridval].Edge_Points[5] = NbPointIsoMap_local;
 2949                     GridVoxelVarPt[(i-1)*maxgrscalemaxgr+(j-1)*maxgridval].NbEdgePoint += 1;
 2950                 }
 2951                 NbPointIsoMap_local++;
 2952             }
 2953         }
 2954 
 2955     /// 6) Case k = nb_depth -1
 2956 
 2957     k = masterthread->XYZgrid -1;
 2958     for(i=0; i < masterthread->XYZgrid; i++)
 2959         for(j=0; j < masterthread->XYZgrid; j++)
 2960         {
 2961             IsoValue_1 = Results[i*maxgrscalemaxgr+j*maxgridval+k];
 2962             // First Case P(i+1)(j)(k)
 2963             if( i != (masterthread->XYZgrid -1) )
 2964             {
 2965                 IsoValue_2 = Results[(i+1)*maxgrscalemaxgr+j*maxgridval+k];
 2966                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 2967                 {
 2968                     // Edge Point computation and  save in IsoPointMap
 2969                     factor = (IsoValue - IsoValue_1)/rapport;
 2970                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i] - factor * masterthread->x_Step[isoindex];
 2971                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j];
 2972                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 2973                     for(int iter=0; iter<7; iter++)
 2974                         NormVertexTabVector.push_back(1.0);
 2975                     NormVertexTabVector.push_back(float(vals[0]));
 2976                     NormVertexTabVector.push_back(float(vals[1]));
 2977                     NormVertexTabVector.push_back(float(vals[2]));
 2978                     // save The reference to this point
 2979                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].Edge_Points[0] = NbPointIsoMap_local;
 2980                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].NbEdgePoint += 1;
 2981                     // The same Point is used in three other Voxels
 2982                     if(j !=0)
 2983                     {
 2984                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k].Edge_Points[4] = NbPointIsoMap_local;
 2985                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k].NbEdgePoint += 1;
 2986                     }
 2987                     if(k !=0)
 2988                     {
 2989                         GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k-1].Edge_Points[2]   = NbPointIsoMap_local;
 2990                         GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k-1].NbEdgePoint += 1;
 2991                     }
 2992                     if(j !=0 && k != 0)
 2993                     {
 2994                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k-1].Edge_Points[6] = NbPointIsoMap_local;
 2995                         GridVoxelVarPt[i*maxgrscalemaxgr+(j-1)*maxgridval+k-1].NbEdgePoint += 1;
 2996                     }
 2997                     NbPointIsoMap_local++;
 2998                 }
 2999             } /// End of if(i != nb_ligne-1)...
 3000 
 3001             // Second Case P(i)(j+1)(k)
 3002             if( j != (masterthread->XYZgrid -1) )
 3003             {
 3004                 IsoValue_2 = Results[i*maxgrscalemaxgr+(j+1)*maxgridval+k];
 3005                 // Edge Point computation and  save in IsoPointMap
 3006                 if(IsoValue_1 * IsoValue_2 <= 0 && (rapport=IsoValue_2 - IsoValue_1) != 0.0)
 3007                 {
 3008                     factor = (IsoValue - IsoValue_1)/rapport;
 3009                     vals[0] = masterthread->xLocal2[isoindex*masterthread->GridVal+i];
 3010                     vals[1] = masterthread->yLocal2[isoindex*masterthread->GridVal+j] - factor * masterthread->y_Step[isoindex];
 3011                     vals[2] = masterthread->zLocal2[isoindex*masterthread->GridVal+k];
 3012                     for(int iter=0; iter<7; iter++)
 3013                         NormVertexTabVector.push_back(1.0);
 3014                     NormVertexTabVector.push_back(float(vals[0]));
 3015                     NormVertexTabVector.push_back(float(vals[1]));
 3016                     NormVertexTabVector.push_back(float(vals[2]));
 3017 
 3018                     // save The reference to this point
 3019                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].Edge_Points[8] = NbPointIsoMap_local;
 3020                     GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k].NbEdgePoint += 1;
 3021 
 3022                     // The same Point is used in three other Voxels
 3023                     if( i !=0 )
 3024                     {
 3025                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval+k].Edge_Points[9]    = NbPointIsoMap_local;
 3026                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval+k].NbEdgePoint += 1;
 3027                     }
 3028                     if( k !=0 )
 3029                     {
 3030                         GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k-1].Edge_Points[11]   = NbPointIsoMap_local;
 3031                         GridVoxelVarPt[i*maxgrscalemaxgr+j*maxgridval+k-1].NbEdgePoint += 1;
 3032                     }
 3033                     if( i !=0 && k !=0 )
 3034                     {
 3035                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval+k-1].Edge_Points[10] = NbPointIsoMap_local;
 3036                         GridVoxelVarPt[(i-1)*maxgrscalemaxgr+j*maxgridval+k-1].NbEdgePoint += 1;
 3037                     }
 3038                     NbPointIsoMap_local++;
 3039                 }
 3040             } /// End of if( j != (nb_colon -1) )...
 3041         }
 3042     NbPointIsoMap =    NbPointIsoMap_local;
 3043     return 1;
 3044 }
 3045 
 3046 ///+++++++++++++++++++++++++++++++++++++++++++++++++++++///
 3047 void Iso3D::SignatureComputation()
 3048 {
 3049     uint I, J, IJK, IPLUSONEJK,
 3050          IJPLUSONEK,  IPLUSONEJPLUSONEK;
 3051     uint maxgrscalemaxgr = masterthread->GridVal*masterthread->GridVal;
 3052     for(uint i=0; i < masterthread->XYZgrid; i++)
 3053     {
 3054         I = i*maxgrscalemaxgr;
 3055         for(uint j=0; j < masterthread->XYZgrid; j++)
 3056         {
 3057             J = I + j*masterthread->GridVal;
 3058             for(uint k=0; k < masterthread->XYZgrid; k++)
 3059             {
 3060                 IJK               = J + k;
 3061                 IPLUSONEJK        = IJK + maxgrscalemaxgr;
 3062                 IJPLUSONEK        = IJK + masterthread->GridVal;
 3063                 IPLUSONEJPLUSONEK = IPLUSONEJK + masterthread->GridVal;
 3064 
 3065                 if(Results[IJK] <= 0) GridVoxelVarPt[IJK].Signature +=1;
 3066 
 3067                 if(i != (masterthread->XYZgrid-1))
 3068                     if(Results[IPLUSONEJK] <= 0) GridVoxelVarPt[IJK].Signature +=2;
 3069 
 3070                 if(i != (masterthread->XYZgrid-1) && k != (masterthread->XYZgrid-1))
 3071                     if(Results[IPLUSONEJK+1] <= 0) GridVoxelVarPt[IJK].Signature +=4;
 3072 
 3073                 if(k != (masterthread->XYZgrid-1))
 3074                     if(Results[IJK+1] <= 0) GridVoxelVarPt[IJK].Signature +=8;
 3075 
 3076                 if(j != (masterthread->XYZgrid-1))
 3077                     if(Results[IJPLUSONEK] <= 0) GridVoxelVarPt[IJK].Signature +=16;
 3078 
 3079                 if(i != (masterthread->XYZgrid-1) && j != (masterthread->XYZgrid-1))
 3080                     if(Results[IPLUSONEJPLUSONEK] <= 0) GridVoxelVarPt[IJK].Signature +=32;
 3081 
 3082                 if(i != (masterthread->XYZgrid-1) && j != (masterthread->XYZgrid-1) && k != (masterthread->XYZgrid-1))
 3083                     if(Results[IPLUSONEJPLUSONEK+1] <= 0) GridVoxelVarPt[IJK].Signature +=64;
 3084 
 3085                 if(j != (masterthread->XYZgrid-1) && k != (masterthread->XYZgrid-1))
 3086                     if(Results[IJPLUSONEK+1] <= 0) GridVoxelVarPt[IJK].Signature +=128;
 3087             }
 3088         }
 3089     }// End if(Grid...
 3090 }