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