rsvs3D  0.0.0
Codes for the c++ implementation of the 3D RSVS
snake.cpp
1 #define _USE_MATH_DEFINES
2 
3 #include "snake.hpp"
4 
5 #include <cmath>
6 #include <iostream>
7 #include <limits>
8 #include <stdexcept>
9 #include <vector>
10 
11 #include "mesh.hpp"
12 #include "meshprocessing.hpp"
13 #include "rsvsutils.hpp"
14 #include "warning.hpp"
15 
16 using namespace std;
17 
18 // Implementation of features
19 
20 // Class Method implementaton
21 
22 // --------------------------------------------------------------------------
23 // Implementatation of snax - snaxedge - snaxsurf
24 // --------------------------------------------------------------------------
25 void snax::disp() const
26 {
27  cout << "snax #" << index << "; d " << d << "; v " << v << "; edgeind " << edgeind << "; fromvert " << fromvert
28  << "; tovert " << tovert << "; isfreeze " << isfreeze << "; orderedge " << orderedge << " | isBorder "
29  << isBorder << endl;
30 }
31 
32 void snaxedge::disp() const
33 {
34  cout << "snaxedge #" << index << " | surfind " << surfind << " "
35  << " | isBorder " << isBorder << " ";
36  normvector.disp();
37 }
38 
39 void snaxsurf::disp() const
40 {
41  cout << "snaxsurf #" << index << index << " | voluind " << voluind << " "
42  << " | isBorder " << isBorder << " ";
43  normvector.disp();
44 }
45 
46 // file i/o
47 void snax::write(FILE *fid) const
48 {
49  fprintf(fid, "%i %lf %lf %i %i %i %i %i \n", index, d, v, edgeind, fromvert, tovert, isfreeze, orderedge);
50 }
51 
52 void snaxedge::write(FILE *fid) const
53 {
54  fprintf(fid, "%i %i %lf %lf %lf \n", index, surfind, normvector(0), normvector(1), normvector(2));
55 }
56 
57 void snaxsurf::write(FILE *fid) const
58 {
59  fprintf(fid, "%i %i %lf %lf %lf \n", index, voluind, normvector(0), normvector(1), normvector(2));
60 }
61 void snax::read(FILE *fid)
62 {
63  fscanf(fid, "%i %lf %lf %i %i %i %i %i ", &index, &d, &v, &edgeind, &fromvert, &tovert, &isfreeze, &orderedge);
64 }
65 
66 void snaxedge::read(FILE *fid)
67 {
68  fscanf(fid, "%i %i %lf %lf %lf ", &index, &surfind, &normvector[0], &normvector[1], &normvector[2]);
69 }
70 
71 void snaxsurf::read(FILE *fid)
72 {
73  fscanf(fid, "%i %i %lf %lf %lf ", &index, &voluind, &normvector[0], &normvector[1], &normvector[2]);
74 }
75 
76 void snaxedge::PrepareForUse()
77 {
78  normvector.PrepareForUse();
79 }
80 
81 void snaxsurf::PrepareForUse()
82 {
83  normvector.PrepareForUse();
84 }
85 
86 #pragma GCC diagnostic push
87 #pragma GCC diagnostic ignored "-Wunused-parameter"
88 void snax::disptree(const mesh &meshin, int n) const
89 {
90  disp();
91  meshin.verts.isearch(index)->disptree(meshin, n);
92 }
93 void snaxedge::disptree(const mesh &meshin, int n) const
94 {
95  disp();
96  meshin.edges.isearch(index)->disptree(meshin, n);
97 }
98 void snaxsurf::disptree(const mesh &meshin, int n) const
99 {
100  disp();
101  meshin.surfs.isearch(index)->disptree(meshin, n);
102 }
103 void snax::disptree(const snake &snakein, int n) const
104 {
105  int i;
106  for (i = 0; i < n; i++)
107  {
108  cout << " ";
109  }
110  disp();
111  for (i = 0; i < n; i++)
112  {
113  cout << " ";
114  }
115  snakein.snakeconn.verts.isearch(index)->disp();
116  if (n > 0)
117  {
118  for (i = 0; unsigned_int(i) < snakein.snakeconn.verts.isearch(index)->edgeind.size(); i++)
119  {
120  snakein.snaxedges.isearch(snakein.snakeconn.verts.isearch(index)->edgeind[i])->disptree(snakein, n - 1);
121  }
122  }
123 }
124 void snaxedge::disptree(const snake &snakein, int n) const
125 {
126  int i;
127  for (i = 0; i < n; i++)
128  {
129  cout << " ";
130  }
131  disp();
132  for (i = 0; i < n; i++)
133  {
134  cout << " ";
135  }
136  snakein.snakeconn.edges.isearch(index)->disp();
137  if (n > 0)
138  {
139  for (i = 0; unsigned_int(i) < snakein.snakeconn.edges.isearch(index)->vertind.size(); i++)
140  {
141  snakein.snaxs.isearch(snakein.snakeconn.edges.isearch(index)->vertind[i])->disptree(snakein, n - 1);
142  }
143  for (i = 0; unsigned_int(i) < snakein.snakeconn.edges.isearch(index)->surfind.size(); i++)
144  {
145  snakein.snaxsurfs.isearch(snakein.snakeconn.edges.isearch(index)->surfind[i])->disptree(snakein, n - 1);
146  }
147  }
148 }
149 void snaxsurf::disptree(const snake &snakein, int n) const
150 {
151  int i;
152  for (i = 0; i < n; i++)
153  {
154  cout << " ";
155  }
156  disp();
157  for (i = 0; i < n; i++)
158  {
159  cout << " ";
160  }
161  snakein.snakeconn.surfs.isearch(index)->disp();
162  if (n > 0)
163  {
164  for (i = 0; unsigned_int(i) < snakein.snakeconn.surfs.isearch(index)->edgeind.size(); i++)
165  {
166  snakein.snaxedges.isearch(snakein.snakeconn.surfs.isearch(index)->edgeind[i])->disptree(snakein, n - 1);
167  }
168  }
169 }
170 
171 void snax::ChangeIndices(int nVert, int nEdge, int nSurf, int nVolu)
172 {
173  index += nVert;
174 }
175 void snaxedge::ChangeIndices(int nVert, int nEdge, int nSurf, int nVolu)
176 {
177  index += nEdge;
178 }
179 void snaxsurf::ChangeIndices(int nVert, int nEdge, int nSurf, int nVolu)
180 {
181  index += nSurf;
182 }
183 void snax::ChangeIndicesSnakeMesh(int nVert, int nEdge, int nSurf, int nVolu)
184 {
185  fromvert += nVert;
186  tovert += nVert;
187  edgeind += nEdge;
188 }
189 void snaxedge::ChangeIndicesSnakeMesh(int nVert, int nEdge, int nSurf, int nVolu)
190 {
191  surfind += nSurf;
192 }
193 void snaxsurf::ChangeIndicesSnakeMesh(int nVert, int nEdge, int nSurf, int nVolu)
194 {
195  voluind += nVolu;
196 }
197 #pragma GCC diagnostic pop
198 
199 // ---------------------------------------------------------------------------
200 // Extensions of SnakStruct for snaxarray
201 // ---------------------------------------------------------------------------
202 
203 bool snaxarray::checkready()
204 {
205  this->readyforuse = SnakStruct<snax>::checkready();
206  this->readyforuse = (this->readyforuse && this->isOrderedOnEdge);
207  return (this->readyforuse);
208 }
209 void snaxarray::PrepareForUse()
210 {
211  SnakStruct<snax>::PrepareForUse();
212  if (isOrderedOnEdge == 0)
213  {
214  this->ReorderOnEdge();
215  }
216  readyforuse = true;
217 }
218 void snaxarray::Concatenate(const snaxarray &other)
219 {
220  SnakStruct<snax>::Concatenate(other);
221  isOrderedOnEdge = 0;
222 }
223 void snaxarray::ForceArrayReady()
224 {
225  SnakStruct<snax>::ForceArrayReady();
226  isOrderedOnEdge = 1;
227 }
228 
229 // -------------------------------------------------------------------------------------------
230 // Implementatation of snake
231 // -------------------------------------------------------------------------------------------
232 
233 void snake::disp() const
234 {
235  snaxs.disp();
236  snaxedges.disp();
237  snaxsurfs.disp();
238  snakeconn.disp();
239  cout << "Snaking Mesh with ";
240  snakemesh()->displight();
241 }
242 void snake::displight() const
243 {
244  cout << "snake: snax " << snaxs.size();
245  cout << "; snaxedges " << snaxedges.size();
246  cout << "; snaxsurfs " << snaxsurfs.size() << endl;
247 
248  cout << "Snake with connectivity ";
249  snakeconn.displight();
250 
251  cout << "Snaking Mesh with ";
252  snakemesh()->displight();
253 }
254 
255 void snake::read(FILE *fid)
256 {
257  snaxs.read(fid);
258  snaxedges.read(fid);
259  snaxsurfs.read(fid);
260 
261  snakeconn.read(fid);
262 }
263 void snake::write(FILE *fid) const
264 {
265  // fprintf(fid,"%i \n",int(borderIsSet));
266  snaxs.write(fid);
267  snaxedges.write(fid);
268  snaxsurfs.write(fid);
269 
270  snakeconn.write(fid);
271 }
272 int snake::read(const char *str)
273 {
274  // Convenience read function taking a single char argument
275  FILE *fid;
276  fid = fopen(str, "r");
277  if (fid != NULL)
278  {
279  this->read(fid);
280  fclose(fid);
281  }
282  else
283  {
284  cout << "File " << str << "Could not be opened to read" << endl;
285  return (1);
286  }
287  return (0);
288 }
289 int snake::write(const char *str) const
290 {
291  FILE *fid;
292  fid = fopen(str, "w");
293  if (fid != NULL)
294  {
295  this->write(fid);
296  fclose(fid);
297  }
298  else
299  {
300  cout << "File " << str << "Could not be opened to write" << endl;
301  return (1);
302  }
303  return (0);
304 }
305 
306 bool snake::isready() const
307 {
308  bool readyforuse = true;
309  readyforuse = readyforuse & snaxs.isready();
310  readyforuse = readyforuse & snaxedges.isready();
311  readyforuse = readyforuse & snaxsurfs.isready();
312  readyforuse = readyforuse & snakeconn.isready();
313  readyforuse = readyforuse & snakemesh()->isready();
314 
315  return (readyforuse);
316 }
317 
318 inline void snake::GetMaxIndex(int *nVert, int *nEdge, int *nSurf, int *nVolu) const
319 {
320  snakeconn.GetMaxIndex(nVert, nEdge, nSurf, nVolu);
321 }
322 
323 void snake::PrepareForUse(bool needOrder)
324 {
325  snaxs.isInMesh = true;
326  snaxedges.isInMesh = true;
327  snaxsurfs.isInMesh = true;
328 
329  snaxs.PrepareForUse();
330  snaxedges.PrepareForUse();
331  snaxsurfs.PrepareForUse();
332  snakeconn.PrepareForUse(needOrder);
333  if (this->snakemesh() != NULL)
334  {
335  snakemesh()->PrepareForUse();
336  }
337  else
338  {
339  cerr << endl
340  << "Warning: Mesh containing snake unset in snake"
341  " object.\n use `SetSnakeMesh(mesh *snakemeshin)` to set attribute"
342  " `snakemesh` to a mesh pointer"
343  << endl;
344  return;
345  }
346 
347  is3D = int(snakemesh()->volus.size()) > 0;
348 
349  if (int(isMeshVertIn.size()) != int(snakemesh()->verts.size()))
350  {
351  isMeshVertIn.assign(snakemesh()->verts.size(), false);
352  }
353 }
354 void snake::HashArray()
355 {
356  // prefer use of PrepareForUse in final code as it includes checks to
357  // avoid unnecessary repetiion of work
358  // Use HashArray when debugging to be sure that fields are correctly hashed
359  snaxs.HashArray();
360  snaxedges.HashArray();
361  snaxsurfs.HashArray();
362  snakeconn.HashArray();
363  snakemesh()->HashArray();
364 }
365 void snake::HashParent()
366 {
367  // prefer use of PrepareForUse in final code as it includes checks to
368  // avoid unnecessary repetiion of work
369  // Use HashArray when debugging to be sure that fields are correctly hashed
370  snaxs.HashParent();
371  snaxedges.HashParent();
372  snaxsurfs.HashParent();
373 }
374 void snake::HashArrayNM()
375 {
376  // prefer use of PrepareForUse in final code as it includes checks to
377  // avoid unnecessary repetiion of work
378  // Use HashArray when debugging to be sure that fields are correctly hashed
379  snaxs.HashArray();
380  snaxedges.HashArray();
381  snaxsurfs.HashArray();
382  snakeconn.HashArray();
383 }
384 void snake::SetMaxIndex()
385 {
386  // prefer use of PrepareForUse in final code as it includes checks to
387  // avoid unnecessary repetiion of work
388 
389  snaxs.SetMaxIndex();
390  snaxedges.SetMaxIndex();
391  snaxsurfs.SetMaxIndex();
392  snakeconn.SetMaxIndex();
393  snakemesh()->SetMaxIndex();
394 }
395 void snake::SetMaxIndexNM()
396 {
397  // prefer use of PrepareForUse in final code as it includes checks to
398  // avoid unnecessary repetiion of work
399 
400  snaxs.SetMaxIndex();
401  snaxedges.SetMaxIndex();
402  snaxsurfs.SetMaxIndex();
403  snakeconn.SetMaxIndex();
404 }
405 void snake::SetLastIndex()
406 {
407  // prefer use of PrepareForUse in final code as it includes checks to
408  // avoid unnecessary repetiion of work
409 
410  snaxs.SetLastIndex();
411  snaxedges.SetLastIndex();
412  snaxsurfs.SetLastIndex();
413  snakeconn.SetLastIndex();
414  snakemesh()->SetLastIndex();
415 }
416 
417 void snake::Init(mesh *snakemeshin, int nSnax, int nEdge, int nSurf, int nVolu)
418 {
419  // Initialise a snake
420  is3D = snakemeshin->volus.size() > 0;
421 
422  snaxs.Init(nSnax);
423  snaxedges.Init(nEdge);
424  snaxsurfs.Init(nSurf);
425 
426  snaxs.isInMesh = true;
427  snaxedges.isInMesh = true;
428  snaxsurfs.isInMesh = true;
429 
430  snakeconn.Init(nSnax, nEdge, nSurf, nVolu);
431 
432  this->SetSnakeMesh(snakemeshin);
433 }
434 
435 void snake::SetSnakeMesh(mesh *snakemeshin)
436 {
437  this->privatesnakemesh = snakemeshin;
438  this->isSetStepLimit = false;
439  this->isMeshVertIn.assign(snakemeshin->verts.size(), false);
440  this->edgeStepLimit.assign(this->snakemesh()->edges.size(), 1.0);
441 }
442 void snake::reserve(int nSnax, int nEdge, int nSurf, int nVolu)
443 {
444  // Initialise a
445 
446  snaxs.reserve(nSnax);
447  snaxedges.reserve(nEdge);
448  snaxsurfs.reserve(nSurf);
449 
450  snaxs.isInMesh = true;
451  snaxedges.isInMesh = true;
452  snaxsurfs.isInMesh = true;
453 
454  snakeconn.reserve(nSnax, nEdge, nSurf, nVolu);
455 }
456 
457 void snake::ChangeIndices(int nVert, int nEdge, int nSurf, int nVolu)
458 {
459  snaxs.ChangeIndices(nVert, nEdge, nSurf, nVolu);
460  snaxedges.ChangeIndices(nVert, nEdge, nSurf, nVolu);
461  snaxsurfs.ChangeIndices(nVert, nEdge, nSurf, nVolu);
462  snakeconn.ChangeIndices(nVert, nEdge, nSurf, nVolu);
463 }
464 void snake::ChangeIndicesSnakeMesh(int nVert, int nEdge, int nSurf, int nVolu)
465 {
466  int ii = 0;
467  for (ii = 0; ii < snaxs.size(); ++ii)
468  {
469  snaxs[ii].ChangeIndicesSnakeMesh(nVert, nEdge, nSurf, nVolu);
470  }
471  for (ii = 0; ii < snaxedges.size(); ++ii)
472  {
473  snaxedges[ii].ChangeIndicesSnakeMesh(nVert, nEdge, nSurf, nVolu);
474  }
475  for (ii = 0; ii < snaxsurfs.size(); ++ii)
476  {
477  snaxsurfs[ii].ChangeIndicesSnakeMesh(nVert, nEdge, nSurf, nVolu);
478  }
479  // Nothing to do for snakeconn
480  // snakeconn.ChangeIndicesSnakeMesh(nVert,nEdge,nSurf,nVolu);
481 
482  snakemesh()->ChangeIndices(nVert, nEdge, nSurf,
483  nVolu); // Maybe this one is a bad idea
484 }
485 
486 int CompareSnakeInternalStatus(const vector<bool> &thisVec, bool thisFlipped, const vector<bool> &otherVec,
487  bool otherFlipped)
488 {
489  /*
490  Checks if a snake is inside another using information from the points inside
491  the snaking mesh.
492  Works by finding counter examples.
493  This is not a complete test and assumes a relatively simple relationship
494  between the two snakes. of either one completely internal or the two side by
495  side.
496  */
497  bool isIndep;
498  int ii, ni;
499  ni = thisVec.size();
500  isIndep = true;
501  for (ii = 0; ii < ni; ++ii)
502  {
503  if ((thisVec[ii] ^ thisFlipped) && (otherVec[ii] ^ otherFlipped))
504  {
505  isIndep = false;
506  break;
507  }
508  }
509 
510  return (int(isIndep));
511 }
512 void snake::Concatenate(const snake &other, int isInternal)
513 {
514  /*
515  isInternal=0 status unknown
516  isInternal=1 other snake is inside this
517  isInternal=2 this snake is inside other
518  isInternal=-1 other and this are independant
519  */
520  int i, ni, isIndep;
521  if (this->snakemesh() == other.snakemesh())
522  {
523  this->snaxs.Concatenate(other.snaxs);
524  this->snaxedges.Concatenate(other.snaxedges);
525  this->snaxsurfs.Concatenate(other.snaxsurfs);
526  this->snakeconn.Concatenate(other.snakeconn);
527  // Maintain the internal status
528  ni = isMeshVertIn.size();
529 
530  if (isInternal == 0)
531  {
532  isIndep = CompareSnakeInternalStatus(isMeshVertIn, ReturnFlip(), other.isMeshVertIn, other.ReturnFlip());
533  }
534  else
535  {
536  isIndep = isInternal <= 0;
537  }
538  // cout << "Indep status " << isIndep << isFlipped << other.isFlipped << " "
539  // ;
540  if (isIndep)
541  {
542  for (i = 0; i < ni; ++i)
543  {
544  isMeshVertIn[i] =
545  ((isMeshVertIn[i] ^ isFlipped) || (other.isMeshVertIn[i] ^ other.isFlipped)) ^ isFlipped;
546  }
547  }
548  else
549  {
550  for (i = 0; i < ni; ++i)
551  {
552  isMeshVertIn[i] =
553  ((isMeshVertIn[i] ^ isFlipped) && (other.isMeshVertIn[i] ^ other.isFlipped)) ^ isFlipped;
554  }
555  }
556 
557  // if(!isFlipped && !other.isFlipped){
558  // for (i=0;i<ni;++i){
559  // isMeshVertIn[i]= isMeshVertIn[i] || other.isMeshVertIn[i];
560  // }
561  // } else if (isFlipped && other.isFlipped) {
562  // for (i=0;i<ni;++i){
563  // isMeshVertIn[i]= isMeshVertIn[i] || other.isMeshVertIn[i];
564  // }
565  // } else if (!isFlipped && other.isFlipped) {
566  // for (i=0;i<ni;++i){
567  // isMeshVertIn[i]= isMeshVertIn[i] ^ other.isMeshVertIn[i];
568  // } //xor
569  // } else if (isFlipped && !other.isFlipped) {
570  // for (i=0;i<ni;++i){
571  // isMeshVertIn[i]= isMeshVertIn[i] ^ other.isMeshVertIn[i];
572  // } //xor
573  // }
574  }
575  else
576  {
577  RSVS3D_ERROR_ARGUMENT("Concatenation of snakes with different base meshes not supported");
578  }
579 }
580 
581 void snake::MakeCompatible_inplace(snake &other) const
582 {
583  // Makes other mesh compatible with this to be
584  // merged without index crashes
585 
586  int nVert, nEdge, nVolu, nSurf, ii;
587 
588  // Define Max indices in current mesh
589  this->GetMaxIndex(&nVert, &nEdge, &nSurf, &nVolu);
590  other.ChangeIndices(nVert, nEdge, nSurf, nVolu);
591  for (ii = 0; ii < other.snaxs.size(); ++ii)
592  {
593  other.snaxs[ii].orderedge = -1;
594  }
595 }
596 snake snake::MakeCompatible(snake other) const
597 {
598  MakeCompatible_inplace(other);
599  return (other);
600 }
601 
602 void snake::VertIsIn(int vertInd, bool isIn)
603 {
604  if (this->snakemesh() != NULL)
605  {
606  int i = this->snakemesh()->verts.find(vertInd);
607  if (i != -1)
608  {
609  isMeshVertIn[i] = isIn;
610  }
611  }
612 }
613 void snake::VertIsIn(vector<int> vertInd, bool isIn)
614 {
615  if (this->snakemesh() != NULL)
616  {
617  int nj = vertInd.size();
618  for (int j = 0; j < nj; j++)
619  {
620  int i = this->snakemesh()->verts.find(vertInd[j]);
621  if (i != -1)
622  {
623  isMeshVertIn[i] = isIn;
624  }
625  }
626  }
627 }
628 
629 void snake::CheckConnectivity() const
630 {
631  int ii, jj, ni, kk, ll, errCount, errTot;
632  vector<int> testSub;
633  bool isErr;
634 
635  testSub.reserve(10);
636 
637  errCount = 0;
638  errTot = 0;
639  ni = int(snaxs.size());
640  for (ii = 0; ii < ni; ++ii)
641  {
642  if (snaxs.countparent(snaxs(ii)->KeyParent()) > 1)
643  {
644  snaxs.findsiblings(snaxs(ii)->KeyParent(), testSub);
645  for (jj = 0; jj < int(testSub.size()); ++jj)
646  {
647  isErr = false;
648  if (snaxedges(ii)->index != snaxedges(testSub[jj])->index)
649  {
650  for (kk = 0; kk < int(snakeconn.verts(ii)->edgeind.size()); kk++)
651  {
652  if ((snakeconn.edges.isearch(snakeconn.verts(ii)->edgeind[kk])->vertind[0] ==
653  snaxs(testSub[jj])->index) ||
654  (snakeconn.edges.isearch(snakeconn.verts(ii)->edgeind[kk])->vertind[1] ==
655  snaxs(testSub[jj])->index))
656  {
657  isErr = true;
658  break;
659  }
660  }
661  }
662  if (isErr)
663  {
664  cerr << endl
665  << " Test Snake Connectivity Error :" << errCount << " snaxel " << snakeconn.verts(ii)->index
666  << " is connected to snaxel " << snaxs(testSub[jj])->index << " list: " << endl;
667  snakeconn.verts(ii)->disp();
668  snakeconn.verts(testSub[jj])->disp();
669  errCount++;
670  }
671  }
672  }
673  }
674  if (errCount > 0)
675  {
676  cerr << "Test Connectivity snax (same edge connectivity) Errors :" << errCount << endl;
677  }
678  errTot += errCount;
679  errCount = 0;
680  ni = int(snaxs.size());
681  /*Snaxel does not have the same number of edges as its parent edge
682  has surface connections*/
683  for (ii = 0; ii < ni; ++ii)
684  {
685  if (snakeconn.verts(ii)->edgeind.size() != snakemesh()->edges.isearch(snaxs(ii)->edgeind)->surfind.size())
686  {
687  cerr << endl
688  << " Test snax number of edges error :" << errCount << " snaxel-vertex " << snakeconn.verts(ii)->index
689  << " is connected to snaxel " << snaxs(ii)->edgeind << " list: " << endl;
690  snakeconn.verts(ii)->disp();
691  snakemesh()->edges.isearch(snaxs(ii)->edgeind)->disp();
692  errCount++;
693  }
694  }
695  if (errCount > 0)
696  {
697  cerr << "Test Connectivity snax (edges in surfs) Errors :" << errCount << endl;
698  }
699  errTot += errCount;
700 
701  errCount = 0;
702  errTot = 0;
703  ni = int(snaxedges.size());
704  for (ii = 0; ii < ni; ++ii)
705  {
706  if (snaxedges.countparent(snaxedges(ii)->KeyParent()) > 1)
707  {
708  snaxedges.findsiblings(snaxedges(ii)->KeyParent(), testSub);
709  for (jj = 0; jj < int(testSub.size()); ++jj)
710  {
711  isErr = false;
712  if (snaxedges(ii)->index != snaxedges(testSub[jj])->index)
713  {
714  for (kk = 0; kk < int(snakeconn.edges(ii)->vertind.size()); kk++)
715  {
716  if ((snakeconn.edges(ii)->vertind[kk] == snakeconn.edges(testSub[jj])->vertind[0]) ||
717  (snakeconn.edges(ii)->vertind[kk] == snakeconn.edges(testSub[jj])->vertind[1]))
718  {
719  isErr = true;
720  break;
721  }
722  }
723  }
724  if (isErr)
725  {
726  cerr << endl
727  << " Test SnakeEdges Connectivity Error :" << errCount << " snaxedge "
728  << snakeconn.edges(ii)->index << " is connected to snaxedge " << snaxedges(testSub[jj])->index
729  << " list: " << endl;
730  snakeconn.edges(ii)->disp();
731  snakeconn.edges(testSub[jj])->disp();
732  errCount++;
733  }
734  }
735  }
736  }
737  if (errCount > 0)
738  {
739  cerr << "Test Connectivity snaxedge (same surf connectivity) Errors :" << errCount << endl;
740  }
741  errTot += errCount;
742 
743  errCount = 0;
744  errTot = 0;
745  ni = int(snaxsurfs.size());
746  for (ii = 0; ii < ni; ++ii)
747  {
748  if (snaxsurfs.countparent(snaxsurfs(ii)->KeyParent()) > 1)
749  {
750  snaxsurfs.findsiblings(snaxsurfs(ii)->KeyParent(), testSub);
751  for (jj = 0; jj < int(testSub.size()); ++jj)
752  {
753  isErr = false;
754  if (snaxsurfs(ii)->index != snaxsurfs(testSub[jj])->index)
755  {
756  for (kk = 0; kk < int(snakeconn.surfs(ii)->edgeind.size()); kk++)
757  {
758  for (ll = 0; ll < int(snakeconn.surfs(testSub[jj])->edgeind.size()); ++ll)
759  {
760  if (snakeconn.surfs(testSub[jj])->edgeind[ll] == snakeconn.surfs(ii)->edgeind[kk])
761  {
762  isErr = true;
763  break;
764  }
765  }
766  if (isErr)
767  {
768  break;
769  }
770  }
771  }
772  if (isErr)
773  {
774  cerr << endl
775  << " Test SnakeSurf Connectivity Error :" << errCount << " snaxsurf "
776  << snakeconn.surfs(ii)->index << " is connected to snaxsurf " << snaxsurfs(testSub[jj])->index
777  << " list: " << endl;
778  snakeconn.surfs(ii)->disp();
779  snakeconn.surfs(testSub[jj])->disp();
780  errCount++;
781  }
782  }
783  }
784  }
785  if (errCount > 0)
786  {
787  cerr << "Test Connectivity snaxsurf (same volu connectivity) Errors :" << errCount << endl;
788  }
789  errTot += errCount;
790 
791  if (errTot > 0)
792  {
793  cerr << errTot << " Total errors were detected in the connectivity list" << endl;
794  }
795 }
796 
797 void snake::OrderEdges()
798 {
799  int nNewSurfs;
800  snaxsurf newSurf;
801  vector<int> newSurfPairs;
802 
803  newSurfPairs = snakeconn.OrderEdges();
804  nNewSurfs = int(newSurfPairs.size() / 2);
805  for (int i = 0; i < nNewSurfs; ++i)
806  {
807  newSurf = *(snaxsurfs.isearch(newSurfPairs[i * 2]));
808  newSurf.index = newSurfPairs[i * 2 + 1];
809  snaxsurfs.push_back(newSurf);
810  }
811  // int nVerts = this->snakeconn.verts.size();
812  // for (int i =0; i<nVerts ; ++i){
813  // this->snakeconn.verts.elems[i].OrderEdges(&(this->snakeconn));
814  // }
815 }
816 
817 // Snake Movement
818 void snake::UpdateDistance(double dt, double maxDstep, bool scaledStep)
819 {
820  int ii;
821  double dStep = 0.0;
822  int nSnax = int(snaxs.size());
823 
824  if (!this->isSetStepLimit && scaledStep)
825  {
826  this->SetEdgeStepLimits();
827  }
828 
829  for (ii = 0; ii < nSnax; ++ii)
830  {
831  double indivMaxDstep = 1.0 * maxDstep;
832  if (scaledStep)
833  {
834  indivMaxDstep = this->SnaxStepLimit(ii) * indivMaxDstep;
835  }
836  if (snaxs[ii].v * dt > indivMaxDstep)
837  {
838  dStep = indivMaxDstep;
839  }
840  else if (snaxs[ii].v * dt < -indivMaxDstep)
841  {
842  dStep = -indivMaxDstep;
843  }
844  else
845  {
846  dStep = snaxs[ii].v * dt;
847  }
848  snaxs[ii].d = snaxs[ii].d + dStep;
849  // Make sure d is in [0.0, 1.0]
850  snaxs[ii].d = min(max(snaxs[ii].d, 0.0), 1.0);
851  }
852 }
853 
854 void snake::UpdateDistance(const vector<double> &dt, double maxDstep, bool scaledStep)
855 {
856  int ii;
857  double dStep = 0.0;
858  if (int(dt.size()) != snaxs.size())
859  {
860  RSVS3D_ERROR_TYPE("Mismatching sizes passed to Update distance of snake"
861  "Snake and Dt had mismatched sizes at snake::UpdateDistance",
862  length_error);
863  }
864 
865  if (!this->isSetStepLimit && scaledStep)
866  {
867  this->SetEdgeStepLimits();
868  }
869 
870  int nSnax = int(snaxs.size());
871  for (ii = 0; ii < nSnax; ++ii)
872  {
873  double indivMaxDstep = 1.0 * maxDstep;
874  if (scaledStep)
875  {
876  indivMaxDstep = this->SnaxStepLimit(ii) * indivMaxDstep;
877  }
878  if (snaxs[ii].v * dt[ii] > indivMaxDstep)
879  {
880  dStep = indivMaxDstep;
881  }
882  else if (snaxs[ii].v * dt[ii] < -indivMaxDstep)
883  {
884  dStep = -indivMaxDstep;
885  }
886  else
887  {
888  dStep = snaxs[ii].v * dt[ii];
889  }
890  snaxs[ii].d = snaxs[ii].d + dStep;
891  // Make sure d is in [0.0, 1.0]
892  snaxs[ii].d = min(max(snaxs[ii].d, 0.0), 1.0);
893  }
894 }
895 
902 void snake::SetEdgeStepLimits()
903 {
904  int nEdges = this->snakemesh()->edges.size();
905  int numParents = this->snakemesh()->meshtree.NumberOfParents();
906  std::vector<int> edgeAccess;
907  HashedVector<int, int> edgeInds;
908  edgeAccess.assign(nEdges, 0);
909  this->edgeStepLimit.assign(nEdges, 0.0);
910  edgeInds.reserve(nEdges);
911 
912  for (int i = 0; i < numParents; ++i)
913  {
914  // go through each Parent volus
915  auto parentmesh = this->snakemesh()->meshtree.ParentPointer(i);
916  // identify all the child volumes
917  // identify the edges they touch
918  int nParentVolus = parentmesh->volus.size();
919  for (int j = 0; j < nParentVolus; ++j)
920  {
921  double minEdgeLength = INFINITY;
922  auto snakeMeshVolus = this->snakemesh()->meshtree.ChildIndices(i, parentmesh->volus(j)->index);
923  edgeInds.clear();
924  for (auto iVolu : snakeMeshVolus)
925  {
926  for (auto iSurf : this->snakemesh()->volus(iVolu)->surfind)
927  {
928  for (auto iEdge : this->snakemesh()->surfs.isearch(iSurf)->edgeind)
929  {
930  if (edgeInds.find(iEdge) == rsvs3d::constants::__notfound)
931  {
932  minEdgeLength = min(minEdgeLength, this->snakemesh()->edges.isearch(iEdge)->GetLength());
933  edgeInds.push_back(iEdge);
934  }
935  }
936  }
937  }
938 
939  for (auto iEdge : edgeInds.vec)
940  {
941  int iEdgePos = this->snakemesh()->edges.find(iEdge);
942  edgeAccess[iEdgePos]++;
943  this->edgeStepLimit[iEdgePos] += minEdgeLength / this->snakemesh()->edges(iEdgePos)->GetLength();
944  }
945  }
946  }
947  for (int i = 0; i < nEdges; ++i)
948  {
949  if (edgeAccess[i] > 1)
950  {
951  this->edgeStepLimit[i] = this->edgeStepLimit[i] / double(edgeAccess[i]);
952  }
953  else if (edgeAccess[i] < 1)
954  {
955  this->edgeStepLimit[i] = 1.0;
956  }
957  }
958  if (numParents > 0)
959  {
960  this->isSetStepLimit = true;
961  }
962 }
963 
964 double snake::SnaxStepLimit(int snaxSub) const
965 {
966  int iEdge = this->snakemesh()->edges.find(this->snaxs(snaxSub)->edgeind);
967 #ifdef SAFE_ACCESS
968  // this->disp();
969  if (iEdge == rsvs3d::constants::__notfound || iEdge >= int(this->edgeStepLimit.size()))
970  {
971  std::cerr << snaxSub << "," << this->snaxs.size() << "," << this->edgeStepLimit.size() << "," << iEdge << ","
972  << this->snaxs(snaxSub)->edgeind;
973  RSVS3D_ERROR("Attemp to access an edge which was not found in "
974  "snake::snakemesh into snake::edgeStepLimit.");
975  }
976 #endif
977  return this->edgeStepLimit[iEdge];
978 }
979 
980 void snake::UpdateCoord()
981 {
982  int fromVertSub, toVertSub;
983  int ii, jj;
984  int nSnax = int(snaxs.size());
985  for (ii = 0; ii < nSnax; ++ii)
986  {
987  fromVertSub = snakemesh()->verts.find(snaxs(ii)->fromvert);
988  toVertSub = snakemesh()->verts.find(snaxs(ii)->tovert);
989 
990  for (jj = 0; jj < 3; ++jj)
991  {
992  snakeconn.verts.elems[ii].coord[jj] =
993  (snakemesh()->verts(toVertSub)->coord[jj] - snakemesh()->verts(fromVertSub)->coord[jj]) * snaxs(ii)->d +
994  snakemesh()->verts(fromVertSub)->coord[jj];
995  }
996  for (auto iEdge : snakeconn.verts(ii)->edgeind)
997  {
998  this->snakeconn.InvalidateEdgeLength(iEdge);
999  }
1000  }
1001 
1002  this->snakeconn.SetEdgeLengths();
1003 }
1004 
1005 void snake::UpdateCoord(const vector<int> &snaxInds)
1006 {
1007  int fromVertSub, toVertSub;
1008  int ii, jj;
1009  for (auto iSnax : snaxInds)
1010  {
1011  ii = this->snaxs.find(iSnax);
1012  fromVertSub = snakemesh()->verts.find(snaxs(ii)->fromvert);
1013  toVertSub = snakemesh()->verts.find(snaxs(ii)->tovert);
1014 
1015  for (jj = 0; jj < 3; ++jj)
1016  {
1017  snakeconn.verts.elems[ii].coord[jj] =
1018  (snakemesh()->verts(toVertSub)->coord[jj] - snakemesh()->verts(fromVertSub)->coord[jj]) * snaxs(ii)->d +
1019  snakemesh()->verts(fromVertSub)->coord[jj];
1020  }
1021  for (auto iEdge : snakeconn.verts(ii)->edgeind)
1022  {
1023  this->snakeconn.InvalidateEdgeLength(iEdge);
1024  }
1025  }
1026 
1027  this->snakeconn.SetEdgeLengths();
1028 }
1029 
1030 // snaxarray extension
1031 
1032 void snaxarray::OrderOnEdge()
1033 {
1034  vector<int> edgeInds, snaxSubs, orderSubs;
1035  vector<double> orderSnax, unorderSnax;
1036  unordered_multimap<double, int> hashOrder;
1037  int isBwd;
1038  int ii, jj, nEdge;
1039 
1040  edgeInds = ConcatenateScalarField(*this, &snax::edgeind, 0, elems.size());
1041  sort(edgeInds);
1042  unique(edgeInds);
1043  for (ii = 0; ii < int(edgeInds.size()); ++ii)
1044  {
1045  nEdge = hashParent.count(edgeInds[ii]);
1046 
1047  if (nEdge > 1)
1048  {
1049  snaxSubs = ReturnDataEqualRange(edgeInds[ii], hashParent);
1050  orderSnax.resize(snaxSubs.size());
1051  for (jj = 0; jj < int(snaxSubs.size()); ++jj)
1052  {
1053  isBwd = elems[snaxSubs[jj]].fromvert > elems[snaxSubs[jj]].tovert;
1054  orderSnax[jj] = isBwd ? (1.0 - elems[snaxSubs[jj]].d) : elems[snaxSubs[jj]].d;
1055  }
1056  unorderSnax = orderSnax;
1057  // This will break with snakes in same place
1058  sort(orderSnax);
1059 
1060  orderSubs = FindSubList(orderSnax, unorderSnax, hashOrder);
1061  for (jj = 0; jj < int(snaxSubs.size()); ++jj)
1062  {
1063  elems[snaxSubs[orderSubs[jj]]].orderedge = jj + 1;
1064  }
1065  }
1066  }
1067 
1068  isOrderedOnEdge = 1;
1069 }
1070 
1071 void snaxarray::ReorderOnEdge()
1072 {
1073  vector<int> edgeInds, snaxSubs;
1074  vector<double> orderSnax;
1075  int isBwd;
1076  int ii, jj, nEdge, maxOrd = 0;
1077  bool needIncrement;
1078 
1079  edgeInds = ConcatenateScalarField(*this, &snax::edgeind, 0, elems.size());
1080  sort(edgeInds);
1081  unique(edgeInds);
1082  for (ii = 0; ii < int(edgeInds.size()); ++ii)
1083  {
1084  nEdge = hashParent.count(edgeInds[ii]);
1085 
1086  if (nEdge > 1)
1087  {
1088  snaxSubs = ReturnDataEqualRange(edgeInds[ii], hashParent);
1089  orderSnax.resize(snaxSubs.size());
1090  needIncrement = false;
1091  maxOrd = 0;
1092  for (jj = 0; jj < int(snaxSubs.size()); ++jj)
1093  {
1094  maxOrd = (maxOrd <= elems[snaxSubs[jj]].orderedge) ? elems[snaxSubs[jj]].orderedge : maxOrd;
1095  }
1096  for (jj = 0; jj < int(snaxSubs.size()); ++jj)
1097  {
1098  if (elems[snaxSubs[jj]].orderedge == -1)
1099  {
1100  isBwd = elems[snaxSubs[jj]].fromvert > elems[snaxSubs[jj]].tovert;
1101 
1102  elems[snaxSubs[jj]].orderedge = (isBwd && (elems[snaxSubs[jj]].d == 0.0))
1103  ? maxOrd + 1
1104  : ((!isBwd && (elems[snaxSubs[jj]].d == 0.0))
1105  ? 0
1106  : ((isBwd && (elems[snaxSubs[jj]].d == 1.0))
1107  ? 0
1108  : ((!isBwd && (elems[snaxSubs[jj]].d == 1.0))
1109  ? maxOrd + 1
1110  : elems[snaxSubs[jj]].orderedge)));
1111  if (elems[snaxSubs[jj]].orderedge == -1)
1112  {
1113  cerr << "Error Order not correctly assigned" << endl;
1114  }
1115  if (elems[snaxSubs[jj]].orderedge == 0)
1116  {
1117  needIncrement = true;
1118  }
1119  }
1120  }
1121  if (needIncrement)
1122  {
1123  for (jj = 0; jj < int(snaxSubs.size()); ++jj)
1124  {
1125  elems[snaxSubs[jj]].orderedge++;
1126  }
1127  }
1128  }
1129  else
1130  {
1131  snaxSubs = ReturnDataEqualRange(edgeInds[ii], hashParent);
1132  elems[snaxSubs[0]].orderedge = 1;
1133  }
1134  }
1135 
1136  isOrderedOnEdge = 1;
1137 }
1138 
1139 void snaxarray::CalculateTimeStepOnEdge(vector<double> &dt, vector<bool> &isSnaxDone, int edgeInd)
1140 {
1141  int nSnax, ii, jj;
1142  vector<int> snaxSubs;
1143  double impactTime;
1144 
1145  snaxSubs = ReturnDataEqualRange(edgeInd, hashParent);
1146  nSnax = snaxSubs.size();
1147  for (ii = 0; ii < nSnax; ++ii)
1148  {
1149  isSnaxDone[snaxSubs[ii]] = true;
1150  for (jj = ii + 1; jj < nSnax; ++jj)
1151  {
1152  impactTime = SnaxImpactDt(elems[snaxSubs[ii]], elems[snaxSubs[jj]]);
1153  if (impactTime >= 0.0)
1154  {
1155  dt[snaxSubs[ii]] = (dt[snaxSubs[ii]] > impactTime) ? impactTime : dt[snaxSubs[ii]];
1156  dt[snaxSubs[jj]] = (dt[snaxSubs[jj]] > impactTime) ? impactTime : dt[snaxSubs[jj]];
1157  }
1158  }
1159  }
1160 }
1161 
1162 double ValidateEquilibrateD(double newD, double oldD, double stepLength)
1163 {
1164  double maxLim = pow(stepLength, 1.5);
1165  double minLim = pow(stepLength, 0.5);
1166  if (oldD < 0.5)
1167  {
1168  if (!isfinite(newD))
1169  {
1170  newD = oldD;
1171  }
1172  else if (newD < minLim)
1173  {
1174  newD = minLim;
1175  }
1176  else if (newD > maxLim)
1177  {
1178  newD = maxLim;
1179  }
1180  }
1181  else
1182  {
1183  if (!isfinite(newD))
1184  {
1185  newD = oldD;
1186  }
1187  else if (newD < (1.0 - maxLim))
1188  {
1189  newD = (1.0 - maxLim);
1190  }
1191  else if (newD > (1.0 - minLim))
1192  {
1193  newD = (1.0 - minLim);
1194  }
1195  }
1196 
1197  return newD;
1198 }
1199 
1200 std::pair<double, double> EquilibrateEdge(snake &snakein, int snaxA, int snaxB, int vertA, int vertB)
1201 {
1202  // calculate centre of edge.
1203  coordvec c, del1, vec;
1204  double dg0;
1205  c = snakein.snakeconn.verts.isearch(snaxA)->coord;
1206  c.add(snakein.snakeconn.verts.isearch(snaxB)->coord);
1207  c.div(2.0);
1208 
1209  del1 = snakein.snakeconn.verts.isearch(vertA)->coord;
1210  c.substract(snakein.snakeconn.verts.isearch(vertB)->coord);
1211  del1.substract(snakein.snakeconn.verts.isearch(vertB)->coord);
1212  vec = del1.cross(c.usedata());
1213 
1214  // c is the plane normal at this point;
1215  double d = vec.dot(snakein.snakeconn.verts.isearch(vertA)->coord);
1216 
1217  // Vertex A
1218  del1 = snakein.snakemesh()->verts.isearch(snakein.snaxs.isearch(snaxA)->tovert)->coord;
1219  del1.substract(snakein.snakemesh()->verts.isearch(snakein.snaxs.isearch(snaxA)->fromvert)->coord);
1220  dg0 = vec.dot(snakein.snakemesh()->verts.isearch(snakein.snaxs.isearch(snaxA)->fromvert)->coord);
1221  double dga = vec.dot(del1.usedata());
1222  double da = (d - dg0) / dga;
1223 
1224  // Vertex B
1225  dg0 = vec.dot(snakein.snakemesh()->verts.isearch(snakein.snaxs.isearch(snaxB)->fromvert)->coord);
1226  del1 = snakein.snakemesh()->verts.isearch(snakein.snaxs.isearch(snaxB)->tovert)->coord;
1227  del1.substract(snakein.snakemesh()->verts.isearch(snakein.snaxs.isearch(snaxB)->fromvert)->coord);
1228  double dgb = vec.dot(del1.usedata());
1229  double db = (d - dg0) / dgb;
1230 
1231  return {da, db};
1232 }
1233 // template <class T=int, int n=2>
1234 void DisplayVectorArray(const std::vector<std::array<int, 2>> &in)
1235 {
1236  int count = in.size();
1237  std::cout << count << ",";
1238  if (!count)
1239  {
1240  return;
1241  }
1242  int count2 = in[0].size();
1243  std::cout << count2 << " : ";
1244  for (int i = 0; i < count; ++i)
1245  {
1246  for (int j = 0; j < count2; ++j)
1247  {
1248  std::cout << in[i][j] << " ";
1249  }
1250  std::cout << " : ";
1251  }
1252  std::cout << " | ";
1253 }
1254 
1255 int SmoothStep_bordersetting(int spawnVert, snake &snakein, double spawnDist)
1256 {
1257  // identify snaxels
1258 
1259  HashedVector<int, int> snaxInds;
1260  double spawnDistTol = spawnDist * 1.1;
1261  snaxInds.reserve(20);
1262  for (auto parentEdge : snakein.snakemesh()->verts.isearch(spawnVert)->edgeind)
1263  {
1264  std::vector<int> snaxIndsTemp;
1265  snaxIndsTemp.reserve(20);
1266  snakein.snaxs.findsiblings(parentEdge, snaxIndsTemp);
1267 
1268  for (auto snaxCandidate : snaxIndsTemp)
1269  {
1270  auto currSnax = snakein.snaxs(snaxCandidate);
1271  if (currSnax->CloseToVertex() == spawnVert &&
1272  snaxInds.find(currSnax->index) == rsvs3d::constants::__notfound &&
1273  (currSnax->d < spawnDistTol || (1 - spawnDistTol) < currSnax->d))
1274  {
1275  snaxInds.push_back(currSnax->index);
1276  }
1277  }
1278  }
1279  if (snaxInds.size() < 2)
1280  {
1281  return 0;
1282  }
1283  HashedVector<int, int> borderSurfs;
1284  for (auto iSnax : snaxInds.vec)
1285  {
1286  for (auto iSurf : snakein.snakeconn.verts.isearch(iSnax)->elmind(snakein.snakeconn, 2))
1287  {
1288  bool isSpawned = false;
1289  for (auto iVert : snakein.snakeconn.surfs.isearch(iSurf)->vertind(snakein.snakeconn))
1290  {
1291  if (iVert == rsvs3d::constants::__notfound)
1292  {
1293  isSpawned = true;
1294  break;
1295  }
1296  }
1297  if (!isSpawned && borderSurfs.find(iSurf) == rsvs3d::constants::__notfound)
1298  {
1299  borderSurfs.push_back(iSurf);
1300  }
1301  }
1302  }
1303 
1304  // Consistancy of this list can be checked making sure all the snaxels come
1305  // from the same vertex or go to the same vertex.
1306  HashedVector<int, int> externalPts;
1307  std::vector<std::array<int, 2>> externalLinks;
1308  externalPts.reserve(snaxInds.size());
1309  externalLinks.reserve(snaxInds.size());
1310  for (auto iSnax : snaxInds.vec)
1311  {
1312  for (auto iEdge : snakein.snakeconn.verts.isearch(iSnax)->edgeind)
1313  {
1314  for (auto iVert : snakein.snakeconn.edges.isearch(iEdge)->vertind)
1315  {
1316  if (snaxInds.find(iVert) == rsvs3d::constants::__notfound)
1317  {
1318  externalPts.push_back(iSnax);
1319  externalLinks.push_back({iEdge, iVert});
1320  }
1321  }
1322  }
1323  }
1324  // Compute the two distances
1325  std::vector<int> edgeProcess;
1326  edgeProcess.reserve(borderSurfs.size());
1327  for (auto border : borderSurfs.vec)
1328  {
1329  bool success = false;
1330  for (auto iEdge : snakein.snakeconn.surfs.isearch(border)->edgeind)
1331  {
1332  auto &vertTemp = snakein.snakeconn.edges.isearch(iEdge)->vertind;
1333  if (snaxInds.find(vertTemp[0]) != rsvs3d::constants::__notfound &&
1334  snaxInds.find(vertTemp[1]) != rsvs3d::constants::__notfound && success)
1335  {
1336  success = false;
1337  edgeProcess.pop_back();
1338  break;
1339  }
1340  else if (snaxInds.find(vertTemp[0]) != rsvs3d::constants::__notfound &&
1341  snaxInds.find(vertTemp[1]) != rsvs3d::constants::__notfound)
1342  {
1343  success = true;
1344  edgeProcess.push_back(iEdge);
1345  }
1346  }
1347  if (!success)
1348  {
1349  edgeProcess.push_back(rsvs3d::constants::__notfound);
1350  // RSVS3D_ERROR_LOGIC("A border edge was not found when looking for it.");
1351  }
1352  }
1353 
1354  std::vector<double> newD;
1355  std::vector<int> countAccess;
1356  newD.assign(snaxInds.size(), 0.0);
1357  countAccess.assign(snaxInds.size(), 0);
1358  int count1 = edgeProcess.size();
1359  for (int i = 0; i < count1; ++i)
1360  {
1361  int actSurf = borderSurfs[i];
1362  int iEdge = edgeProcess[i];
1363  if (iEdge == rsvs3d::constants::__notfound)
1364  {
1365  continue;
1366  }
1367  // pass the edge index, the two snaxels, the two vertices
1368  int snaxA = snakein.snakeconn.edges.isearch(iEdge)->vertind[0];
1369  int snaxB = snakein.snakeconn.edges.isearch(iEdge)->vertind[1];
1370 
1371  int posA = rsvs3d::constants::__notfound;
1372  int posB = rsvs3d::constants::__notfound;
1373  bool success = false;
1374 
1375  for (auto iA : externalPts.findall(snaxA))
1376  {
1377  // Find the one where they all share the same surface
1378  posA = iA;
1379  auto testEdge = snakein.snakeconn.edges.isearch(externalLinks[posA][0]);
1380  if (actSurf == testEdge->surfind[0] || actSurf == testEdge->surfind[1])
1381  {
1382  success = true;
1383  break;
1384  }
1385  }
1386  if (!success)
1387  {
1388  std::cerr << snaxA << " " << posA << " | ";
1389  DisplayVector(externalPts.vec);
1390  DisplayVector(borderSurfs.vec);
1391  DisplayVector(edgeProcess);
1392  DisplayVectorArray(externalLinks);
1393  std::cerr << std::endl;
1394  snakein.snakeconn.edges.isearch(iEdge)->disp();
1395  snakein.snakeconn.surfs.isearch(actSurf)->disp();
1396  snakein.snakeconn.verts.isearch(snaxA)->disp();
1397  RSVS3D_ERROR_LOGIC("Pos A could not be found");
1398  }
1399  success = false;
1400  for (auto iA : externalPts.findall(snaxB))
1401  {
1402  // Find the one where they all share the same surface
1403  posB = iA;
1404  if (borderSurfs.find(actSurf) != rsvs3d::constants::__notfound)
1405  {
1406  auto testEdge = snakein.snakeconn.edges.isearch(externalLinks[posA][0]);
1407  if (actSurf == testEdge->surfind[0] || actSurf == testEdge->surfind[1])
1408  {
1409  success = true;
1410  break;
1411  }
1412  }
1413  }
1414  if (!success)
1415  {
1416  RSVS3D_ERROR_LOGIC("Pos B could not be found");
1417  }
1418 
1419  int vertA = externalLinks[posA][1];
1420  int vertB = externalLinks[posB][1];
1421  auto tempD = EquilibrateEdge(snakein, snaxA, snaxB, vertA, vertB);
1422  posA = snaxInds.find(snaxA);
1423  posB = snaxInds.find(snaxB);
1424  ValidateEquilibrateD(tempD.first, snakein.snaxs.isearch(snaxA)->d, spawnDist);
1425  // newD[posA] += tempD.first;
1426  newD[posA] += ValidateEquilibrateD(tempD.first, snakein.snaxs.isearch(snaxA)->d, spawnDist);
1427  countAccess[posA]++;
1428  // newD[posB] += tempD.second;
1429  newD[posB] += ValidateEquilibrateD(tempD.second, snakein.snaxs.isearch(snaxB)->d, spawnDist);
1430  countAccess[posB]++;
1431  }
1432  // How to set all the values not set by the border algorithm?
1433 
1434  // Set the values
1435  int count = snaxInds.size();
1436  int noissue = 0;
1437  for (int i = 0; i < count; ++i)
1438  {
1439  if (countAccess[i] > 0)
1440  {
1441  if ((newD[i] < 0.0 || countAccess[i] != 2 || newD[i] / countAccess[i] > 1.0) && false)
1442  {
1443  DisplayVector(newD);
1444  DisplayVector(countAccess);
1445  std::cout << endl;
1446  DisplayVector(externalPts.vec);
1447  DisplayVector(snaxInds.vec);
1448  DisplayVector(borderSurfs.vec);
1449  DisplayVector(edgeProcess);
1450 
1451  RSVS3D_ERROR_NOTHROW("Negative d was calculated");
1452  }
1453  else
1454  {
1455  noissue++;
1456  snakein.snaxs[snakein.snaxs.find(snaxInds[i], true)].d = newD[i] / double(countAccess[i]);
1457  }
1458  }
1459  }
1460  return (noissue);
1461 }
1462 
1463 namespace rsvs3d
1464 {
1465 namespace smoothsnaking
1466 {
1467 
1477 void FindSnaxInds(HashedVector<int, int> &snaxInds, int spawnVert, const snake &snakein, double spawnDist)
1478 {
1479  double spawnDistTol = spawnDist * 1.1;
1480  snaxInds.reserve(20);
1481  for (auto parentEdge : snakein.snakemesh()->verts.isearch(spawnVert)->edgeind)
1482  {
1483  std::vector<int> snaxIndsTemp;
1484  snaxIndsTemp.reserve(20);
1485  snakein.snaxs.findsiblings(parentEdge, snaxIndsTemp);
1486 
1487  for (auto snaxCandidate : snaxIndsTemp)
1488  {
1489  auto currSnax = snakein.snaxs(snaxCandidate);
1490  if (currSnax->CloseToVertex() == spawnVert &&
1491  snaxInds.find(currSnax->index) == rsvs3d::constants::__notfound &&
1492  (currSnax->d < spawnDistTol || (1.0 - spawnDistTol) < currSnax->d))
1493  {
1494  snaxInds.push_back(currSnax->index);
1495  }
1496  }
1497  }
1498 }
1499 
1500 int LaplacianSmoothingDirection(const snake &snakein, const vert *vertsmooth, coordvec &smoothDir)
1501 {
1502  return VertexLaplacianVector(snakein.snakeconn, vertsmooth, smoothDir, true);
1503 }
1504 } // namespace smoothsnaking
1505 } // namespace rsvs3d
1506 
1524 int SmoothStep_sphere(int spawnVert, snake &snakein, double spawnDist)
1525 {
1526  // identify snaxels
1527 
1528  HashedVector<int, int> snaxInds;
1529  rsvs3d::smoothsnaking::FindSnaxInds(snaxInds, spawnVert, snakein, spawnDist);
1530  if (snaxInds.size() < 2)
1531  {
1532  return 0;
1533  }
1534 
1535  // Consistancy of this list can be checked making sure all the snaxels come
1536  // from the same vertex or go to the same vertex.
1537 
1538  std::vector<int> externalEdges, externalPts, externalSurfsPairs, externalEdgesOrder, internalPts,
1539  internalPtsOrdered;
1540  externalPts.reserve(snaxInds.size());
1541  externalEdges.reserve(snaxInds.size());
1542  internalPts.reserve(snaxInds.size());
1543  externalSurfsPairs.reserve(snaxInds.size() * 2);
1544  int k = 0;
1545  for (auto iSnax : snaxInds.vec)
1546  {
1547  for (auto iEdge : snakein.snakeconn.verts.isearch(iSnax)->edgeind)
1548  {
1549  auto currEdge = snakein.snakeconn.edges.isearch(iEdge);
1550  for (auto iVert : currEdge->vertind)
1551  {
1552  if (snaxInds.find(iVert) == rsvs3d::constants::__notfound)
1553  {
1554  internalPts.push_back(iSnax);
1555  externalEdges.push_back(iEdge);
1556  externalPts.push_back(iVert);
1557  externalSurfsPairs.push_back(currEdge->surfind[0]);
1558  externalSurfsPairs.push_back(currEdge->surfind[1]);
1559  externalEdgesOrder.push_back(k++);
1560  }
1561  }
1562  }
1563  }
1564  int retValue = OrderList(externalEdgesOrder, externalSurfsPairs, false, true);
1565  int nVerts = externalEdgesOrder.size();
1566  if (retValue != rsvs3d::constants::ordering::ordered)
1567  {
1568  RSVS3D_ERROR_NOTHROW("Failed to order the list of edges.");
1569  return 0;
1570  }
1571  if (nVerts < 2)
1572  {
1573  RSVS3D_ERROR_LOGIC("Only one point was left, cannot define direction.");
1574  }
1575  std::vector<const std::vector<double> *> vecPts;
1576  vecPts.reserve(nVerts);
1577  internalPtsOrdered.reserve(nVerts);
1578  for (int i = 0; i < nVerts; ++i)
1579  {
1580  vecPts.push_back(&(snakein.snakeconn.verts.isearch(externalPts[externalEdgesOrder[i]])->coord));
1581  }
1582  auto &centre = snakein.snakemesh()->verts.isearch(spawnVert)->coord;
1583  std::tuple<coordvec, double> centreData;
1584  try
1585  {
1586  centreData = VertexNormal(centre, vecPts);
1587  }
1588  catch (exception const &ex)
1589  {
1590  // This case can happen and just needs to be handled it is the result of
1591  // all the points being the same. An alternative could be to compute the
1592  // normal using the directrion of the snaxels.
1593  return 0;
1594  }
1595  coordvec &vertNormal = get<0>(centreData);
1596  if (rsvs3d::utils::IsAproxEqual(vertNormal.GetNorm(), 0.0))
1597  {
1598  // RSVS3D_ERROR("Vertex normal could not be computed: returned 0 norm.");
1599  return 0;
1600  }
1601  // Figure out if the normal is inward facing or outward facing.
1602  // Order the internal edges
1603  internalPtsOrdered.reserve(nVerts);
1604  int v1 = rsvs3d::constants::__notfound;
1605  int v2 = rsvs3d::constants::__notfound;
1606  for (int i = 0; i < nVerts; ++i)
1607  {
1608  internalPtsOrdered.push_back(internalPts[externalEdgesOrder[i]]);
1609  }
1610  v1 = internalPts[externalEdgesOrder[0]];
1611  int edgeAfterV2 = rsvs3d::constants::__notfound;
1612  for (int i = 1; i < nVerts; ++i)
1613  {
1614  if (v1 != internalPtsOrdered[i])
1615  {
1616  v2 = internalPtsOrdered[i];
1617  edgeAfterV2 = externalEdges[externalEdgesOrder[i]];
1618  break;
1619  }
1620  }
1621  if (v2 == rsvs3d::constants::__notfound)
1622  {
1623  RSVS3D_ERROR_LOGIC("Only one point was found, cannot define direction.");
1624  }
1625  // Extract a surface that is between two points (on an edge)
1626  int edgeInd = snakein.snakeconn.EdgeFromVerts(v1, v2);
1627  if (edgeInd == rsvs3d::constants::__notfound)
1628  {
1629  RSVS3D_ERROR_LOGIC("The two points do not share an edge.");
1630  }
1631  // Compare the list orders
1632  unique(internalPtsOrdered);
1633  if (internalPtsOrdered.back() == internalPtsOrdered[0])
1634  {
1635  internalPtsOrdered.pop_back();
1636  }
1637  if (internalPtsOrdered.size() < 2)
1638  {
1639  RSVS3D_ERROR_LOGIC("Only one point was left, cannot define direction.");
1640  }
1641 
1642  int surfInd = snakein.snakeconn.SurfFromEdges(edgeInd, edgeAfterV2);
1643 #ifdef RSVS_DIAGNOSTIC
1644  std::cerr << std::endl;
1645  std::cerr << "snaxInds ";
1646  DisplayVector(snaxInds.vec);
1647  std::cerr << std::endl;
1648  std::cerr << "internalPtsOrdered ";
1649  DisplayVector(internalPtsOrdered);
1650  std::cerr << std::endl;
1651  std::cerr << "centre ";
1652  DisplayVector(centre);
1653  std::cerr << std::endl;
1654  std::cerr << "outer surfPts ";
1655  DisplayVector(internalPts);
1656  std::cerr << std::endl;
1657  std::cerr << "vertNormal ";
1658  DisplayVector(vertNormal.usedata());
1659  std::cerr << std::endl;
1660  snakein.snakeconn.surfs.isearch(surfInd)->disp();
1661  snakein.snakeconn.edges.isearch(edgeInd)->disp();
1662 #endif
1663 
1664  snakein.snakeconn.surfs.isearch(surfInd)->OrderedVerts(&snakein.snakeconn, internalPts);
1665  int flipMultiplier =
1666  meshhelp::NormalShouldFlip(internalPts, v1, v2, snakein.snakeconn.surfs.isearch(surfInd)->voluind, false);
1667  if (flipMultiplier == -1)
1668  {
1669  vertNormal.flipsign();
1670  }
1671  vertNormal.Normalize(); // Should be the outfacing vertex normal
1672 
1673 #ifdef RSVS_DIAGNOSTIC
1674  std::cerr << " isSameOrder " << isSameOrder << " flipMultiplier " << flipMultiplier << std::endl;
1675 
1676 #endif
1694  int radiusStep = 1;
1695 
1696  if (!snakein.isMeshVertIn[snakein.snakemesh()->verts.find(spawnVert)])
1697  {
1698  // radiusStep = -1;
1699  vertNormal.flipsign();
1700  }
1701  double maxOffset = 0.0;
1702  for (auto iEdge : externalEdges)
1703  {
1704  double maxOffsetTemp = snakein.snakeconn.edges.isearch(iEdge)->Length(snakein.snakeconn);
1705  maxOffset = max(maxOffsetTemp, maxOffset);
1706  }
1707  double minEdgeLength = 0.0;
1708  for (auto iSnax : snaxInds.vec)
1709  {
1710  int iEdge = snakein.snaxs.isearch(iSnax)->edgeind;
1711  double minTemp = snakein.snakemesh()->edges.isearch(iEdge)->Length(*snakein.snakemesh());
1712  minEdgeLength = max(minTemp, minEdgeLength);
1713  }
1714 
1715  double offsetMultiplier = 1.0 - (2.0 * get<1>(centreData) / (2.0 * M_PI));
1716  if (offsetMultiplier < 0.0)
1717  {
1718  return 0;
1719  }
1720  offsetMultiplier = min(max(offsetMultiplier, 0.0), 1.0);
1721  double offset = offsetMultiplier * maxOffset;
1722  vertNormal.mult(offset);
1723  // vertNormal.add(centre); // This is now the sphere centre
1730  // std::array<double, 2> IntersectLineSphere(const coordvec &lineVec,
1731  // vertNormal, sphereRadius)
1732  double sphereRadius = offset + radiusStep * minEdgeLength * spawnDist;
1733  std::vector<double> newD;
1734  coordvec lineVec;
1735  newD.reserve(snaxInds.vec.size());
1736  double maxNew = 0.0;
1737  double infVal = -INFINITY;
1738  bool allNegErr = false;
1739  auto errValFunc = [&](coordvec &lineVec) -> double { return spawnDist * minEdgeLength / lineVec.GetNorm(); };
1740  for (auto iSnax : snaxInds.vec)
1741  {
1742  lineVec = centre;
1743  int farVert = snakein.snaxs.isearch(iSnax)->CloseToVertex(true);
1744  lineVec.substractfrom(snakein.snakemesh()->verts.isearch(farVert)->coord);
1745  auto ds = IntersectLineSphere(lineVec, vertNormal, sphereRadius);
1746  if (radiusStep == -1 && rsvs3d::sign(ds[0]) == 1 && rsvs3d::sign(ds[1]) == 1)
1747  {
1748  newD.push_back(min(ds[0], ds[1]));
1749  }
1750  else if (!isfinite(ds[0]) && ds[0] < -1e+300)
1751  {
1752  newD.push_back(infVal);
1753  }
1754  else if ((rsvs3d::sign(ds[0]) == rsvs3d::sign(ds[1])) || !isfinite(ds[0]))
1755  {
1756  std::cerr << "ds : " << ds[0] << " " << ds[1] << std::endl;
1757  if (!isfinite(ds[0]))
1758  {
1759  std::cerr << std::endl;
1760  int count = vecPts.size();
1761  DisplayVector(centre);
1762  for (int i = 0; i < count; ++i)
1763  {
1764  std::cerr << std::endl;
1765  DisplayVector(*vecPts[i]);
1766  }
1767  RSVS3D_ERROR("The intersection returned Nans.");
1768  }
1769  else if (rsvs3d::sign(ds[0]) == 1)
1770  {
1771  RSVS3D_ERROR_NOTHROW("The intersection returned only positive numbers.");
1772  }
1773  else if (rsvs3d::sign(ds[0]) == -1)
1774  {
1775  RSVS3D_ERROR_NOTHROW("The intersection returned only negative numbers.");
1776  allNegErr = true;
1777  }
1778  newD.push_back(errValFunc(lineVec));
1779  }
1780  else
1781  {
1782  int dUse = (rsvs3d::sign(ds[1]) + 1) / 2;
1783  if (fabs(ds[dUse]) <= fabs(ds[1 - dUse]))
1784  {
1785  newD.push_back(ds[dUse]);
1786  }
1787  else
1788  {
1789  // if this case is triggered it is a sign that the layout
1790  // of snaxels does not match the surface and that the resulting
1791  // surface is not going to be flat.
1792  newD.push_back(errValFunc(lineVec));
1793  }
1794  }
1795 #ifdef RSVS_DIAGNOSTIC
1796  std::cerr << "ds : " << ds[0] << " " << ds[1] << std::endl;
1797 #endif
1798 
1799  maxNew = max(maxNew, newD.back());
1800  }
1801 #ifdef RSVS_DIAGNOSTIC
1802  std::cerr << radiusStep << std::endl << "newD : ";
1803  DisplayVector(newD);
1804  std::cerr << std::endl;
1805 #endif
1806  if (allNegErr)
1807  {
1808  return 0;
1809  }
1810 
1811  double scaleInfinite = 10.0;
1812  double scaleVals = spawnDist / (maxNew);
1813  double minD = spawnDist / (scaleInfinite);
1814  double maxD = scaleInfinite * spawnDist;
1815  int noissue = 0;
1816  int count = snaxInds.vec.size();
1817 
1818  auto edgeScale = [&](coordvec &lineVec) -> double { return minEdgeLength / lineVec.GetNorm(); };
1819  for (int i = 0; i < count; ++i)
1820  {
1821  lineVec = centre;
1822  int farVert = snakein.snaxs.isearch(snaxInds[i])->CloseToVertex(true);
1823  lineVec.substractfrom(snakein.snakemesh()->verts.isearch(farVert)->coord);
1824  double edgeScaleCurr = edgeScale(lineVec);
1825 
1826  newD[i] = scaleVals * newD[i];
1827  newD[i] = (isfinite(newD[i]) ? newD[i] : maxD * edgeScaleCurr);
1828  // newD[i] = (isfinite(newD[i]) && newD[i]<=maxD ? newD[i] :
1829  // maxD*edgeScaleCurr);
1830  newD[i] = max(minD * edgeScaleCurr, newD[i]);
1831  }
1832  // DisplayVector(newD);
1833  for (int i = 0; i < count; ++i)
1834  {
1835  int j = snakein.snaxs.find(snaxInds[i], true);
1836  noissue++;
1837  snakein.snaxs[j].d = snakein.snaxs(j)->fromvert == spawnVert ? newD[i] : 1.0 - newD[i];
1838  snakein.snaxs[j].ValidateDistance(snakein);
1839  }
1840 
1841  return (noissue);
1842 }
1843 
1844 int SmoothStep_planematching(int spawnVert, snake &snakein, double spawnDist)
1845 {
1846  // identify snaxels
1847 
1848  HashedVector<int, int> snaxInds;
1849  double spawnDistTol = spawnDist * 1.1;
1850  snaxInds.reserve(20);
1851  for (auto parentEdge : snakein.snakemesh()->verts.isearch(spawnVert)->edgeind)
1852  {
1853  std::vector<int> snaxIndsTemp;
1854  snaxIndsTemp.reserve(20);
1855  snakein.snaxs.findsiblings(parentEdge, snaxIndsTemp);
1856 
1857  for (auto snaxCandidate : snaxIndsTemp)
1858  {
1859  auto currSnax = snakein.snaxs(snaxCandidate);
1860  if (currSnax->CloseToVertex() == spawnVert &&
1861  snaxInds.find(currSnax->index) == rsvs3d::constants::__notfound &&
1862  (currSnax->d < spawnDistTol || (1 - spawnDistTol) < currSnax->d))
1863  {
1864  snaxInds.push_back(currSnax->index);
1865  }
1866  }
1867  }
1868  if (snaxInds.size() < 2)
1869  {
1870  return 0;
1871  }
1872  // Consistancy of this list can be checked making sure all the snaxels come
1873  // from the same vertex or go to the same vertex.
1874  HashedVector<int, int> externalPts;
1875  std::vector<std::array<int, 2>> externalLinks;
1876  externalPts.reserve(snaxInds.size());
1877  externalLinks.reserve(snaxInds.size());
1878  for (auto iSnax : snaxInds.vec)
1879  {
1880  for (auto iEdge : snakein.snakeconn.verts.isearch(iSnax)->edgeind)
1881  {
1882  for (auto iVert : snakein.snakeconn.edges.isearch(iEdge)->vertind)
1883  {
1884  if (snaxInds.find(iVert) == rsvs3d::constants::__notfound)
1885  {
1886  externalPts.push_back(iSnax);
1887  externalLinks.push_back({iEdge, iVert});
1888  }
1889  }
1890  }
1891  }
1892  // Compute the two distances
1893  std::vector<int> edgeProcess;
1894  edgeProcess.reserve(externalPts.size());
1895  for (auto iSnax : snaxInds.vec)
1896  {
1897  for (auto iEdge : snakein.snakeconn.verts.isearch(iSnax)->edgeind)
1898  {
1899  int iVert1 = snakein.snakeconn.edges.isearch(iEdge)->vertind[0];
1900  int iVert2 = snakein.snakeconn.edges.isearch(iEdge)->vertind[1];
1901  if (externalPts.find(iVert1) != rsvs3d::constants::__notfound &&
1902  externalPts.find(iVert2) != rsvs3d::constants::__notfound)
1903  {
1904  edgeProcess.push_back(iEdge);
1905  }
1906  }
1907  }
1908 
1909  std::vector<double> newD;
1910  newD.assign(externalPts.size(), 0.0);
1911 
1912  for (auto iEdge : edgeProcess)
1913  {
1914  // pass the edge index, the two snaxels, the two vertices
1915  int snaxA = snakein.snakeconn.edges.isearch(iEdge)->vertind[0];
1916  int snaxB = snakein.snakeconn.edges.isearch(iEdge)->vertind[1];
1917  int vertA = externalLinks[externalPts.find(snaxA)][1];
1918  int vertB = externalLinks[externalPts.find(snaxB)][1];
1919  auto tempD = EquilibrateEdge(snakein, snaxA, snaxB, vertA, vertB);
1920  newD[externalPts.find(snaxA)] += tempD.first;
1921  newD[externalPts.find(snaxB)] += tempD.second;
1922  }
1923  // How to set all the values not set by the border algorithm?
1924 
1925  // Set the values
1926  int noissue = 0;
1927  int count = externalPts.size();
1928  for (int i = 0; i < count; ++i)
1929  {
1930  if (newD[i] < __DBL_EPSILON__)
1931  {
1932  DisplayVector(newD);
1933  DisplayVector(externalPts.vec);
1934  DisplayVector(edgeProcess);
1935 
1936  RSVS3D_ERROR_NOTHROW("Negative d was calculated");
1937  }
1938  snakein.snaxs[snakein.snaxs.find(externalPts[i], true)].d = newD[i] / 2.0;
1939  noissue++;
1940  }
1941  return noissue;
1942 }
1943 
1944 int SmoothStep_itersmooth(int spawnVert, snake &snakein, double spawnDist)
1945 {
1946  // Steps:
1947  // + Define the distance bounds that might be performed
1948  // + Gather the snaxels that will be smoothed
1949  // + Implement cotangent laplacian smoothing
1950  int nIter = 5;
1951 
1952  double stepDamping = 2.0 / double(nIter + 1);
1953 
1954  HashedVector<int, int> snaxInds;
1955  snaxInds.clear();
1956  rsvs3d::smoothsnaking::FindSnaxInds(snaxInds, spawnVert, snakein, spawnDist);
1957  if (snaxInds.size() < 2)
1958  {
1959  return 0;
1960  }
1961 
1962  // Order the edges around each snaxel that needs to be analysed.
1963  bool succesOrdering = true;
1964  for (auto iVert : snaxInds.vec)
1965  {
1966  int vertOrdered = snakein.snakeconn.OrderVertexEdges(iVert);
1967  succesOrdering = (vertOrdered == rsvs3d::constants::ordering::ordered);
1968  }
1969  if (!succesOrdering)
1970  { // if ordering failed stop the smoothing process.
1971  RSVS3D_ERROR_NOTHROW("Failed to order vertex");
1972  return 0;
1973  }
1974  //
1975  coordvec smoothDir, snaxDir;
1976  std::vector<double> dSteps;
1977  int nSnax = snaxInds.vec.size();
1978  dSteps.reserve(nSnax);
1979 
1980  for (int stepNum = 0; stepNum < nIter; ++stepNum)
1981  {
1982  dSteps.clear();
1983  double maxD = -INFINITY;
1984  double minD = INFINITY;
1985  double sumD = 0.0;
1986  double minEdgeLength = INFINITY;
1987  double maxEdgeLength = -INFINITY;
1988  double sumEdgeL = 0.0;
1989  for (auto iVert : snaxInds.vec)
1990  {
1991  rsvs3d::smoothsnaking::LaplacianSmoothingDirection(snakein, snakein.snakeconn.verts.isearch(iVert),
1992  smoothDir);
1993  snakein.snakemesh()->VerticesVector(spawnVert, snakein.snaxs.isearch(iVert)->CloseToVertex(true), snaxDir);
1994  minEdgeLength = min(minEdgeLength, snaxDir.GetNorm());
1995  maxEdgeLength = max(maxEdgeLength, snaxDir.GetNorm());
1996  sumEdgeL += snaxDir.GetNorm();
1997  double edgeL = snaxDir.Normalize();
1998  double dStep = snaxDir.dot(smoothDir.usedata()) / edgeL;
1999  dStep = isfinite(dStep) ? dStep : 0.0;
2000  dSteps.push_back(dStep);
2001 
2002  sumD += dStep;
2003  maxD = max(dStep, maxD);
2004  minD = min(dStep, minD);
2005  }
2006  double maxAbs = max(max(fabs(maxD), fabs(minD)), spawnDist);
2007  // double edgeRatio = minEdgeLength/maxEdgeLength;
2008  double stepMultiplier = (spawnDist * stepDamping) / maxAbs;
2009 
2010  for (int i = 0; i < nSnax; ++i)
2011  {
2012  int j = snakein.snaxs.find(snaxInds[i], true);
2013  double l = snakein.edgeStepLimit[snakein.snakemesh()->edges.find(snakein.snaxs(j)->edgeind)];
2014  double tempStep = (snakein.snaxs(j)->fromvert == spawnVert ? dSteps[i] : -dSteps[i]) * stepMultiplier * l;
2015  double prevD = snakein.snaxs[j].d;
2016  snakein.snaxs[j].d += tempStep;
2017  stringstream errstr;
2018  errstr << "Attempted to take a step of length " << tempStep << " snaxel->d: " << prevD;
2019  errstr << " Step " << stepNum << " of " << nIter;
2020  errstr << " maxD " << maxD << " minD " << minD << " l " << l << "dSteps : ";
2021  PrintVector(dSteps, errstr);
2022  errstr << std::endl;
2023  try
2024  {
2025  snakein.snaxs[j].ValidateDistance(snakein);
2026  }
2027  catch (const exception &exc)
2028  {
2029  RSVS3D_ERROR_NOTHROW(errstr.str().c_str());
2030  snakein.snaxs[j].d =
2031  snakein.snaxs[j].d < 0.5 ? spawnDist / double(nIter) : 1.0 - spawnDist / double(nIter);
2032  }
2033 #ifdef RSVS_DIAGNOSTIC_RESOLVED
2034  std::cout << errstr.str();
2035 #endif
2036  }
2037  snakein.snaxs.PrepareForUse();
2038  snakein.UpdateCoord(snaxInds.vec);
2039  // std::cout << meanD << "," << stepMultiplier << ","
2040  // << maxAbs << " | ";
2041  }
2042  double totD = 0.0;
2043  for (auto iSnax : snaxInds.vec)
2044  {
2045  totD += snakein.snaxs.isearch(iSnax)->d < 0.5 ? snakein.snaxs.isearch(iSnax)->d
2046  : 1.0 - snakein.snaxs.isearch(iSnax)->d;
2047  }
2048  std::cout << " totD " << totD / snaxInds.size();
2049  return snaxInds.size();
2050 }
2051 
2052 int SmoothStep(int spawnVert, snake &snakein, double spawnDist, std::string str)
2053 {
2054  static bool errHit = false;
2055 
2056  if (rsvs3d::utils::IsAproxEqual(spawnDist, 0.0))
2057  {
2058  return 0;
2059  }
2060  if (str.compare("none") == 0 || str.compare("") == 0)
2061  {
2062  return 0;
2063  }
2064  else if (str.compare("laplacian") == 0)
2065  {
2066  return SmoothStep_itersmooth(spawnVert, snakein, spawnDist);
2067  }
2068  else if (str.compare("planar") == 0)
2069  {
2070  return SmoothStep_planematching(spawnVert, snakein, spawnDist);
2071  }
2072  else if (str.compare("sphere") == 0)
2073  {
2074  return SmoothStep_sphere(spawnVert, snakein, spawnDist);
2075  }
2076  else
2077  {
2078  if (!errHit)
2079  {
2080  errHit = true;
2081  stringstream errstr;
2082  errstr << "unrecognised smooth snaking step type: '" << str << "' defaulting to 'none'.";
2083  RSVS3D_ERROR_NOTHROW(errstr.str().c_str());
2084  }
2085  }
2086  return 0;
2087 }
2088 
2089 void snake::TakeSpawnStep(int minIndex, double stepLength)
2090 {
2091  int nSnax = this->snaxs.size();
2092  std::vector<int> spawnVerts;
2093  spawnVerts.reserve(nSnax);
2094  for (int i = 0; i < nSnax; ++i)
2095  {
2096  if (this->snaxs(i)->index > minIndex)
2097  {
2098  this->snaxs.elems[i].TakeSpawnStep(*this, stepLength);
2099  }
2100  }
2101 }
2102 
2103 void snake::TakeSmoothSpawnStep(int minIndex, double stepLength, std::string smoothStep)
2104 {
2105  int nSnax = this->snaxs.size();
2106  std::vector<int> spawnVerts;
2107  spawnVerts.reserve(nSnax);
2108  for (int i = 0; i < nSnax; ++i)
2109  {
2110  if (this->snaxs(i)->index > minIndex)
2111  {
2112  spawnVerts.push_back(this->snaxs(i)->CloseToVertex());
2113  }
2114  }
2115  sort(spawnVerts);
2116  unique(spawnVerts);
2117  int numSuccess = 0;
2118  for (auto spawnVert : spawnVerts)
2119  {
2120  numSuccess += SmoothStep(spawnVert, *this, stepLength, smoothStep);
2121  this->snaxs.ForceArrayReady();
2122  }
2123 #ifdef RSVS_DEBUG
2124  std::cout << std::endl << "Number of successful smoothed steps : " << numSuccess << std::endl;
2125 #endif
2126 }
2127 
2128 void snax::TakeSpawnStep(snake &snakein, double stepLength)
2129 {
2130  if (rsvs3d::utils::IsAproxEqual(stepLength, 0.0))
2131  {
2132  return;
2133  }
2134  if (stepLength < -__DBL_EPSILON__ || stepLength >= 0.5)
2135  {
2136  RSVS3D_ERROR_ARGUMENT("The spawn step must be 0.5>=x>=0");
2137  }
2138  int closeVert = this->fromvert;
2139  if (this->d > 0.5)
2140  { // if spawned in reverse
2141  closeVert = this->tovert;
2142  stepLength = -stepLength;
2143  }
2144 
2145  double minEdgeLength = INFINITY;
2146  for (auto edgeInd : snakein.snakemesh()->verts.isearch(closeVert)->edgeind)
2147  {
2148  double testLength = snakein.snakemesh()->edges.isearch(edgeInd)->Length(*snakein.snakemesh());
2149  minEdgeLength = testLength < minEdgeLength ? testLength : minEdgeLength;
2150  }
2151  stepLength =
2152  stepLength * minEdgeLength / snakein.snakemesh()->edges.isearch(this->edgeind)->Length(*snakein.snakemesh());
2153  int nEdge = snakein.snaxs.countparent(this->edgeind);
2154  this->d += stepLength;
2155  if (this->d < 0 || this->d > 1.0)
2156  {
2157  RSVS3D_ERROR_LOGIC("Distance should be between 0 and 1");
2158  }
2159  if (nEdge == 1)
2160  {
2161  return;
2162  }
2163  else if (nEdge > 1)
2164  {
2165  int nChanges = 0;
2166  std::vector<int> snaxSubs;
2167  snakein.snaxs.findsiblings(this->edgeind, snaxSubs);
2168  for (auto snaxtest : snaxSubs)
2169  {
2170  auto tempSnax = snakein.snaxs(snaxtest);
2171  if (this->index != tempSnax->index)
2172  {
2173  if (!((tempSnax->d < 0.5) ^ (this->d < 0.5)))
2174  {
2175  if (fabs(tempSnax->d - 0.5) > fabs(this->d))
2176  {
2177  this->d = tempSnax->d;
2178  nChanges++;
2179  }
2180  }
2181  }
2182  }
2183  if (nChanges > 0)
2184  {
2185  std::cout << " " << nChanges << "," << this->d << ",";
2186  }
2187  }
2188  else
2189  {
2190  RSVS3D_ERROR_RANGE("Parent hashing returned a number <=0");
2191  }
2192 }
2193 
2194 void snax::ValidateDistance(snake &snakein)
2195 {
2196  int nEdge = snakein.snaxs.countparent(this->edgeind);
2197  if (this->d < 0 || this->d > 1.0)
2198  {
2199  stringstream errstr;
2200  errstr << "Distance should be between 0 and 1 was " << this->d;
2201  RSVS3D_ERROR_LOGIC(errstr.str().c_str());
2202  }
2203  if (nEdge == 1)
2204  {
2205  return;
2206  }
2207  else if (nEdge > 1)
2208  {
2209  static std::vector<int> snaxSubs;
2210  snaxSubs.clear();
2211  snakein.snaxs.findsiblings(this->edgeind, snaxSubs);
2212  for (auto snaxtest : snaxSubs)
2213  {
2214  auto tempSnax = snakein.snaxs(snaxtest);
2215  if (this->index != tempSnax->index)
2216  {
2217  // Indicates if snakes are moving in different directions along the edge
2218  // Impacts how to test equality between distances
2219  bool snaxDifferentDir = this->fromvert != tempSnax->fromvert;
2220  // Checks if the edge is being travelled on in reversed or direct
2221  // order.
2222  bool inOrder = tempSnax->fromvert < tempSnax->tovert;
2223  // Checks what the ordering is along that edge
2224  bool tempIsAfter = this->orderedge < tempSnax->orderedge;
2225  // Those three booleans provide all the necessary information to
2226  // know if the snaxels have crossed over and what operations to apply
2227  if (snaxDifferentDir)
2228  {
2229  this->d = 1.0 - this->d;
2230  }
2231  // (inOrder ^ tempIsAfter) -> tells us to test superiority or
2232  // inferiority
2233  // 1 1 0 -> this > temp (modify cases)
2234  // 1 0 1 -> this < temp
2235  // 0 1 1 -> this < temp
2236  // 0 0 0 -> this > temp
2237  if ((inOrder ^ tempIsAfter) && this->d < tempSnax->d)
2238  {
2239  this->d = tempSnax->d;
2240  }
2241  else if (!(inOrder ^ tempIsAfter) && this->d > tempSnax->d)
2242  {
2243  this->d = tempSnax->d;
2244  }
2245 
2246  if (snaxDifferentDir)
2247  {
2248  this->d = 1.0 - this->d;
2249  }
2250  }
2251  }
2252  }
2253  else
2254  {
2255  RSVS3D_ERROR_RANGE("Parent hashing returned a number <=0");
2256  }
2257 }
2258 
2259 void snax::Direction(const snake &snakein, coordvec &dir) const
2260 {
2261  snakein.snakemesh()->VerticesVector(this->fromvert, this->tovert, dir);
2262 }
2263 
2264 // Snake operations
2265 
2266 void snake::ForceCloseContainers()
2267 {
2268  this->snakeconn.ForceCloseContainers();
2269  if (!Check3D())
2270  {
2271  snaxsurfs.elems.clear();
2272  snaxsurfs.Init(snakeconn.surfs.size());
2273  snaxsurfs.PopulateIndices();
2274  snaxsurfs.HashArray();
2275  this->SetSnaxSurfs();
2276  }
2277 }
2278 
2279 void snake::CalculateTimeStep(vector<double> &dt, double dtDefault, double distDefault)
2280 {
2281  int nSnax, ii, nEdge;
2282  vector<bool> isSnaxDone;
2283  double newMindDt = dtDefault;
2284  nSnax = snaxs.size();
2285  isSnaxDone.assign(nSnax, false);
2286 
2287  if (!snaxs.checkready())
2288  {
2289  snaxs.PrepareForUse();
2290  }
2291  if (this->snaxDistanceLimit_conserveShape)
2292  {
2293  for (int i = 0; i < nSnax; ++i)
2294  {
2295  if (fabs(snaxs(i)->v > __DBL_EPSILON__))
2296  {
2297  newMindDt = distDefault / fabs(snaxs(i)->v);
2298  dtDefault = newMindDt < dtDefault ? newMindDt : dtDefault;
2299  }
2300  }
2301  }
2302 
2303  dt.assign(nSnax, dtDefault);
2304  dt.resize(nSnax);
2305 
2306  for (ii = 0; ii < nSnax; ++ii)
2307  {
2308  if (!isSnaxDone[ii])
2309  {
2310  nEdge = snaxs.countparent(snaxs(ii)->edgeind);
2311  if (nEdge == 1)
2312  {
2313  isSnaxDone[ii] = true;
2314  }
2315  else if (nEdge > 1)
2316  {
2317  snaxs.CalculateTimeStepOnEdge(dt, isSnaxDone, snaxs(ii)->edgeind);
2318  }
2319  else
2320  {
2321  cerr << "Error: hashParent was not up to date" << endl;
2322  cerr << "Error in " << __PRETTY_FUNCTION__ << endl;
2323  RSVS3D_ERROR_RANGE("Incorrect hash table provided");
2324  }
2325  }
2326  }
2327 }
2328 
2329 void snake::Flip()
2330 {
2331  int temp;
2332  for (int ii = 0; ii < snaxs.size(); ii++)
2333  {
2334  temp = snaxs[ii].fromvert;
2335  snaxs[ii].fromvert = snaxs(ii)->tovert;
2336  snaxs[ii].tovert = temp;
2337 
2338  snaxs[ii].d = 1.0 - snaxs[ii].d;
2339  snaxs[ii].v = -snaxs[ii].v;
2340  }
2341  isFlipped = !isFlipped;
2342  snaxs.ForceArrayReady();
2343 }
2344 
2345 double SnaxImpactDt(const snax &snax1, const snax &snax2)
2346 {
2347  int isSameDir;
2348  double dD, dV;
2349 
2350  isSameDir = snax1.fromvert == snax2.fromvert;
2351 
2352  dD = ((1.0 * !isSameDir) + (1.0 + (-2.0) * !isSameDir) * snax2.d) - snax1.d;
2353  dV = (1.0 + (-2.0) * !isSameDir) * snax2.v - snax1.v;
2354 
2355  if (rsvs3d::utils::IsAproxEqual(dD, 0.0))
2356  {
2357  return (0.0);
2358  }
2359  if (rsvs3d::utils::IsAproxEqual(dV, 0.0))
2360  {
2361  return (-1.0);
2362  }
2363 
2364  return (-dD / dV);
2365 }
2366 
2367 void snake::SnaxImpactDetection(vector<int> &isImpact)
2368 {
2369  // TODO make an approximate arrival
2370  int nSnax, ii, nEdge;
2371  vector<bool> isSnaxDone;
2372 
2373  nSnax = snaxs.size();
2374  isSnaxDone.assign(nSnax, false);
2375  isImpact.clear();
2376  isImpact.reserve(nSnax * 2);
2377 
2378  if (!snaxs.checkready())
2379  {
2380  snaxs.PrepareForUse();
2381  }
2382 
2383  for (ii = 0; ii < nSnax; ++ii)
2384  {
2385  if (!isSnaxDone[ii])
2386  {
2387  nEdge = snaxs.countparent(snaxs(ii)->edgeind);
2388  if (nEdge == 1)
2389  {
2390  if (rsvs3d::utils::IsAproxEqual(snaxs(ii)->d, 0.0) && (snaxs(ii)->v <= 0.0))
2391  {
2392  isImpact.push_back(snaxs(ii)->index);
2393  isImpact.push_back(-1);
2394  }
2395  else if (rsvs3d::utils::IsAproxEqual(snaxs(ii)->d, 1.0) && (snaxs(ii)->v >= 0.0))
2396  {
2397  isImpact.push_back(snaxs(ii)->index);
2398  isImpact.push_back(-2);
2399  }
2400  isSnaxDone[ii] = true;
2401  }
2402  else if (nEdge > 1)
2403  {
2404  snaxs.DetectImpactOnEdge(isImpact, isSnaxDone, snaxs(ii)->edgeind);
2405  }
2406  else
2407  {
2408  cerr << "Error: hashParent was not up to date" << endl;
2409  cerr << "Error in " << __PRETTY_FUNCTION__ << endl;
2410  RSVS3D_ERROR_RANGE("Incorrect hash table provided");
2411  }
2412  }
2413  }
2414 }
2415 void snake::SnaxAlmostImpactDetection(vector<int> &isImpact, double dDlim)
2416 {
2417  int i;
2418  int nVerts, nSnax, vertInd;
2419  vector<int> vertSnaxCloseCnt, vertSnaxMinInd;
2420  vector<double> vertSnaxMinD;
2421  double d;
2422  bool flag;
2423 
2424  nVerts = snakemesh()->verts.size();
2425  nSnax = snaxs.size();
2426 
2427  // isImpact.clear();
2428  // isImpact.reserve(nSnax*2);
2429 
2430  vertSnaxCloseCnt.assign(nVerts, 0);
2431  vertSnaxMinInd.assign(nVerts, 0);
2432  vertSnaxMinD.assign(nVerts, 1.0);
2433 
2434  /*Find the closest snaxel to each vertex and the number of close
2435  snaxels*/
2436  for (i = 0; i < nSnax; ++i)
2437  {
2438  flag = false;
2439  if ((snaxs(i)->d < dDlim) && (snaxs(i)->v <= 0.0) && !rsvs3d::utils::IsAproxEqual(snaxs(i)->d, 0.0))
2440  {
2441  flag = true;
2442  vertInd = snaxs(i)->fromvert;
2443  d = snaxs(i)->d;
2444  }
2445  else if (((1.0 - snaxs(i)->d) < dDlim) && (snaxs(i)->v >= 0.0) &&
2446  !rsvs3d::utils::IsAproxEqual(snaxs(i)->d, 1.0))
2447  {
2448  flag = true;
2449  vertInd = snaxs(i)->tovert;
2450  d = 1 - snaxs(i)->d;
2451  }
2452  if (flag)
2453  {
2454  vertInd = snakemesh()->verts.find(vertInd);
2455  vertSnaxCloseCnt[vertInd]++;
2456 
2457  if (vertSnaxMinD[vertInd] >= d)
2458  {
2459  vertSnaxMinD[vertInd] = d;
2460  vertSnaxMinInd[vertInd] = i;
2461  }
2462  }
2463  }
2464  // TODO: Need to add a check for the order of the snaxel on the edge.
2465  // if it is not 1 or the largest this action will cause crossovers.
2466  for (i = 0; i < nVerts; ++i)
2467  {
2468  if (vertSnaxCloseCnt[i] >= 2)
2469  {
2470  isImpact.push_back(snaxs(vertSnaxMinInd[i])->index);
2471  if (snaxs(vertSnaxMinInd[i])->d < dDlim)
2472  {
2473  // Collapse on the from vertex
2474  snaxs[vertSnaxMinInd[i]].d = 0.0;
2475  isImpact.push_back(-1);
2476  }
2477  else
2478  {
2479  // Collapse on the to vertex
2480  snaxs[vertSnaxMinInd[i]].d = 1.0;
2481  isImpact.push_back(-2);
2482  }
2483  }
2484  }
2485  snaxs.PrepareForUse();
2486 }
2487 
2488 void snaxarray::DetectImpactOnEdge(vector<int> &isImpact, vector<bool> &isSnaxDone, int edgeInd)
2489 {
2490  int nSnax, ii, jj, dOrd;
2491  vector<int> snaxSubs;
2492  double impactTime;
2493  snaxSubs = ReturnDataEqualRange(edgeInd, hashParent);
2494  nSnax = snaxSubs.size();
2495  for (ii = 0; ii < nSnax; ++ii)
2496  {
2497  isSnaxDone[snaxSubs[ii]] = true;
2498 
2499  for (jj = ii + 1; jj < nSnax; ++jj)
2500  {
2501  impactTime = SnaxImpactDt(elems[snaxSubs[ii]], elems[snaxSubs[jj]]);
2502  dOrd = abs(elems[snaxSubs[ii]].orderedge - elems[snaxSubs[jj]].orderedge);
2503  if (dOrd == 1 && rsvs3d::utils::IsAproxEqual(impactTime, 0.0))
2504  {
2505  isImpact.push_back((elems[snaxSubs[ii]].index));
2506  isImpact.push_back((elems[snaxSubs[jj]].index));
2507 
2508  isImpact.push_back((elems[snaxSubs[jj]].index));
2509  isImpact.push_back((elems[snaxSubs[ii]].index));
2510  }
2511  }
2512 
2513  if (rsvs3d::utils::IsAproxEqual(elems[snaxSubs[ii]].d, 0.0) && (elems[snaxSubs[ii]].v <= 0.0) &&
2514  elems[snaxSubs[ii]].orderedge == 1)
2515  {
2516  isImpact.push_back(elems[snaxSubs[ii]].index);
2517  isImpact.push_back(-1);
2518  }
2519  else if (rsvs3d::utils::IsAproxEqual(elems[snaxSubs[ii]].d, 1.0) && (elems[snaxSubs[ii]].v >= 0.0) &&
2520  elems[snaxSubs[ii]].orderedge == nSnax)
2521  {
2522  isImpact.push_back(elems[snaxSubs[ii]].index);
2523  isImpact.push_back(-2);
2524  }
2525  }
2526 }
2527 
2528 void snake::OrientFaces()
2529 {
2530  /*
2531  Orients either surfaces or edges depending.
2532  */
2533  if (this->Check3D())
2534  {
2535  this->OrientSurfaceVolume();
2536  }
2537  else
2538  {
2539  this->OrientEdgeSurface();
2540  }
2541 }
2542 void snake::OrientEdgeSurface()
2543 {
2544  cerr << "Warning: not coded yet in " << __PRETTY_FUNCTION__ << endl;
2545 }
2546 void snake::OrientSurfaceVolume()
2547 {
2548  // Orders the surf.voluind [c0 c1] such that the surface normal vector points
2549  // from c0 to c1
2550  // This is done by using the surface normals and checking they go towards
2551  // the centre of the cell
2552 
2553  int nBlocks, ii, jj, ni, nj, kk, ll;
2554  vector<int> surfOrient;
2555  vector<bool> isFlip;
2556  double dotProd;
2557  coordvec centreVolu, normalVec;
2558 
2559  nBlocks = snakeconn.OrientRelativeSurfaceVolume(surfOrient);
2560  isFlip.assign(nBlocks, false);
2561  //========================================
2562  // Select direction using coordinate geometry
2563  // use a surface
2564 
2565  for (ii = 1; ii <= nBlocks; ii++)
2566  {
2567  jj = -1;
2568  nj = surfOrient.size();
2569  do
2570  {
2571  jj++;
2572  while (jj < nj && ii != abs(surfOrient[jj]))
2573  {
2574  jj++;
2575  }
2576  if (jj == nj)
2577  { // if the orientation cannot be defined
2578  dotProd = 1.0;
2579  kk = 0;
2580  // cerr << "Warning: Cell orientations could not be computed " << endl;
2581  // cerr << " in " << __PRETTY_FUNCTION__ << endl;
2582  break;
2583  }
2584  kk = snakeconn.surfs(jj)->voluind[0] == 0;
2585  ll = snaxs.isearch(snakeconn.edges.isearch(snakeconn.surfs(jj)->edgeind[0])->vertind[0])->tovert;
2586  centreVolu = snakemesh()->verts.isearch(ll)->coord;
2587  ll = snaxs.isearch(snakeconn.edges.isearch(snakeconn.surfs(jj)->edgeind[0])->vertind[0])->fromvert;
2588  centreVolu.substract(snakemesh()->verts.isearch(ll)->coord);
2589 
2590  normalVec = snakeconn.CalcPseudoNormalSurf(snakeconn.surfs(jj)->index);
2591 
2592  dotProd = centreVolu.dot(normalVec.usedata());
2593  } while (!isfinite(dotProd) || (fabs(dotProd) < numeric_limits<double>::epsilon()));
2594 
2595  isFlip[ii - 1] = (((dotProd < 0.0) && (kk == 0)) || ((dotProd > 0.0) && (kk == 1)));
2596  }
2597  if (nBlocks > 0)
2598  {
2599  ni = surfOrient.size();
2600  for (ii = 0; ii < ni; ++ii)
2601  {
2602  if (isFlip[abs(surfOrient[ii]) - 1])
2603  {
2604  snakeconn.surfs.elems[ii].FlipVolus();
2605  }
2606  }
2607  }
2608 }
2609 void snake::AssignInternalVerts()
2610 {
2611  /*
2612  Assigns internal vertices to the snake based on the snakemesh() connectivity.
2613  */
2614  int ii, ni;
2615  vector<int> vertBlock;
2616 
2617  FindBlockSnakeMeshVerts(vertBlock);
2618  this->isFlipped = false;
2619  ni = vertBlock.size();
2620  for (ii = 0; ii < ni; ++ii)
2621  {
2622  isMeshVertIn[ii] = vertBlock[ii] > 0;
2623  }
2624 }
2625 int snake::FindBlockSnakeMeshVerts(vector<int> &vertBlock) const
2626 {
2627  // int mesh::ConnectedVertex(vector<int> &vertBlock) const{
2628  // Fills a vector with a number for each vertex corresponding to a
2629  // group of connected edges it is part of , can be used close surfaces in 2D
2630  // or volumes in 3D. Uses a flood fill with queue method
2631 
2632  int nVertExplored, nVerts, nBlocks, nCurr, nEdgesCurr, ii, jj, kk, nSnaxs, ll, nl;
2633  vector<bool> vertStatus, snaxStatus; // 1 explored 0 not explored
2634 
2635  vector<int> currQueue, nextQueue,
2636  tempSnax; // Current and next queues of indices
2637 
2638  // Preparation of the arrays;
2639  nVerts = snakemesh()->verts.size();
2640  nSnaxs = snaxs.size();
2641  nBlocks = 0;
2642  nVertExplored = 0;
2643 
2644  snaxStatus.assign(nSnaxs, false);
2645  vertStatus.assign(nVerts, false);
2646  vertBlock.assign(nVerts, 0);
2647  currQueue.reserve(nVerts / 2);
2648  nextQueue.reserve(nVerts / 2);
2649 
2650  // While Loop, while not all vertices explored
2651  while (nVertExplored < nSnaxs)
2652  {
2653  // if currQueue is empty start new block
2654  if (currQueue.size() < 1)
2655  {
2656  // cout << "Block " << nBlocks << " - " << nVertExplored << " - " <<
2657  // nVerts << endl;
2658  ii = 0;
2659  while (ii < nSnaxs && snaxStatus[ii])
2660  {
2661  ii++;
2662  }
2663  if (snaxStatus[ii])
2664  {
2665  cerr << "Error starting point for loop not found despite max number of "
2666  "vertex not reached"
2667  << endl;
2668  cerr << "Error in " << __PRETTY_FUNCTION__ << endl;
2669  RSVS3D_ERROR_RANGE(" : Starting point for block not found");
2670  }
2671  currQueue.push_back(snakemesh()->verts.find(snaxs(ii)->fromvert));
2672  nBlocks++;
2673  }
2674  // Explore current queue
2675  nCurr = currQueue.size();
2676  // cout << "nCurr " << nCurr << " nVertExplored " << nVertExplored <<
2677  // " nSnax " << nSnaxs << " " << nBlocks << endl;
2678  for (ii = 0; ii < nCurr; ++ii)
2679  {
2680  if (!vertStatus[currQueue[ii]])
2681  {
2682  vertBlock[currQueue[ii]] = nBlocks;
2683  nEdgesCurr = snakemesh()->verts(currQueue[ii])->edgeind.size();
2684  for (jj = 0; jj < nEdgesCurr; ++jj)
2685  {
2686  // only follow edge if no snaxel exists on it
2687  if (snaxs.countparent(snakemesh()->verts(currQueue[ii])->edgeind[jj]) < 1)
2688  {
2689  kk = int(
2690  snakemesh()->edges.isearch(snakemesh()->verts(currQueue[ii])->edgeind[jj])->vertind[0] ==
2691  snakemesh()->verts(currQueue[ii])->index);
2692  nextQueue.push_back(snakemesh()->verts.find(
2693  snakemesh()->edges.isearch(snakemesh()->verts(currQueue[ii])->edgeind[jj])->vertind[kk]));
2694 #ifdef SAFE_ALGO
2695  if (snakemesh()->verts.find(snakemesh()
2696  ->edges.isearch(snakemesh()->verts(currQueue[ii])->edgeind[jj])
2697  ->vertind[kk]) == -1)
2698  {
2699  cerr << "Edge index: " << snakemesh()->verts(currQueue[ii])->edgeind[jj] << " vertex index:"
2700  << snakemesh()
2701  ->edges.isearch(snakemesh()->verts(currQueue[ii])->edgeind[jj])
2702  ->vertind[kk]
2703  << endl;
2704  cerr << "Edge connected to non existant vertex" << endl;
2705  cerr << "Error in " << __PRETTY_FUNCTION__ << endl;
2706  RSVS3D_ERROR_RANGE(" : Vertex not found");
2707  }
2708 #endif
2709  }
2710  else
2711  {
2712  snaxs.findsiblings(snakemesh()->verts(currQueue[ii])->edgeind[jj], tempSnax);
2713  nl = tempSnax.size();
2714  for (ll = 0; ll < nl; ++ll)
2715  {
2716  if (snakemesh()->verts(currQueue[ii])->index == snaxs(tempSnax[ll])->fromvert)
2717  {
2718  nVertExplored += int(!snaxStatus[tempSnax[ll]]);
2719  snaxStatus[tempSnax[ll]] = true;
2720  }
2721  }
2722  }
2723  }
2724  vertStatus[currQueue[ii]] = true;
2725  }
2726  }
2727 
2728  // Reset current queue and set to next queue
2729  currQueue.clear();
2730  currQueue.swap(nextQueue);
2731  }
2732  return (nBlocks);
2733 }
2734 grid::limits snake::Scale(const grid::limits &newSize)
2735 {
2736  auto boundBox = this->snakemesh()->BoundingBox();
2737  // Scales the snakemesh() to the correct size
2738  auto transformation = this->snakemesh()->Scale(newSize);
2739  // Uses the same transformation on the snakeconn
2740  this->snakeconn.LinearTransform(transformation);
2741  // Applies the transformation to the family of the snakemesh()
2742  // ie: the VOS mesh.
2743  this->snakemesh()->LinearTransformFamily(transformation);
2744 
2745  return boundBox;
2746 }
2747 
2748 std::vector<double> snake::MoveDirections() const
2749 {
2750  std::vector<double> coords;
2751  int nVert = this->snaxs.size();
2752  coords.reserve(nVert * 3);
2753  coordvec lapVec;
2754 
2755  grid::coordlist neighCoord;
2756 
2757  for (int i = 0; i < nVert; ++i)
2758  {
2759  this->snaxs(i)->Direction(*this, lapVec);
2760  for (int j = 0; j < 3; ++j)
2761  {
2762  coords.push_back(lapVec(j));
2763  }
2764  }
2765  return coords;
2766 }
2767 
2768 // -------------------------------------------------------------------------------------------
2769 // TEST CODE
2770 // -------------------------------------------------------------------------------------------
2771 
2772 // in snakstruct_test.cpp
coordvec
Handles the use and norm of a vector for which the norm and the unit value might be needed.
Definition: mesh.hpp:114
mesh
Class for mesh handling.
Definition: mesh.hpp:592
mesh::SurfFromEdges
int SurfFromEdges(int e1, int e2, int repetitionBehaviour=-1) const
Returns the index of the surface connecting two edges.
Definition: mesh.cpp:5872
snake
Definition: snake.hpp:83
snake::SetEdgeStepLimits
void SetEdgeStepLimits()
Sets the relative edge step limits.
Definition: snake.cpp:902
snax
Definition: snake.hpp:192
vert
Class for a vertex in a mesh.
Definition: mesh.hpp:448
mesh.hpp
Provides all the mesh tools used for the generation of 3D grids and geometries.
grid::coordlist
std::vector< const std::vector< double > * > coordlist
Defines a list of coordinates.
Definition: mesh.hpp:87
VertexNormal
std::tuple< coordvec, double > VertexNormal(const std::vector< double > &centre, const grid::coordlist &vecPts)
Calculates the vertex normal weighted by surface angle partitions.
Definition: mesh.cpp:747
grid::transformation
std::array< std::array< double, 3 >, 3 > transformation
Defines a linear transformation to the mesh where for each dimension: {new minimum,...
Definition: mesh.hpp:85
OrderList
int OrderList(std::vector< int > &edgeind, const std::vector< int > &edge2Vert, bool warn=true, bool errout=true, const std::vector< int > *edgeIndOrigPtr=NULL)
Orders a list of elements defined by pairs of indices.
Definition: mesh.cpp:3449
meshhelp::NormalShouldFlip
int NormalShouldFlip(const std::vector< int > orderedList, int elm1, int elm2, const std::vector< int > &voluind, bool innerComparison)
Answers if a normal should be flipped to point towards the outside of a volume.
Definition: mesh.cpp:777
meshprocessing.hpp
Tools for the mathematical processing of meshes.
IntersectLineSphere
std::array< double, 2 > IntersectLineSphere(const coordvec &lineVec, const coordvec &offset, double sphereRadius)
Calculates the parametric position along a line at which it intersects a sphere.
Definition: meshprocessing.cpp:1070
rsvs3d
Namespace for general purpose tools of the RSVS project.
Definition: snake.cpp:1464
rsvs3d::sign
int sign(T val)
Returns the sign of a type comparable to 0.
Definition: warning.hpp:215
rsvsutils.hpp
Provides various utility functions for the RSVS operation.
snake.hpp
Provides the core restricted surface snake container.
warning.hpp
Provides the error and warning system used by the RSVS3D project.
RSVS3D_ERROR_LOGIC
#define RSVS3D_ERROR_LOGIC(M)
Throw a logic_error.
Definition: warning.hpp:139
RSVS3D_ERROR_NOTHROW
#define RSVS3D_ERROR_NOTHROW(M)
Generic rsvs warning.
Definition: warning.hpp:120
RSVS3D_ERROR
#define RSVS3D_ERROR(M)
Throw generic rsvs errors.
Definition: warning.hpp:113
RSVS3D_ERROR_RANGE
#define RSVS3D_ERROR_RANGE(M)
Throw a range_error.
Definition: warning.hpp:173
RSVS3D_ERROR_ARGUMENT
#define RSVS3D_ERROR_ARGUMENT(M)
Throw a invalid_argument.
Definition: warning.hpp:148
RSVS3D_ERROR_TYPE
#define RSVS3D_ERROR_TYPE(M, T)
Throw a specific error type.
Definition: warning.hpp:130