UnityMol  0.9.6-875
UnityMol viewer / In developement
MarchingCubesRec.cs
Go to the documentation of this file.
1 
66 using UnityEngine;
67 using System.Collections.Generic;
68 using System;
69 using Molecule.Model;
70 
71 public class MarchingCubesRec {
72 
73  private int ncellsX; // dimensions of the grid
74  private int ncellsY;
75  private int ncellsZ;
76 // private float toleranceValue; // range around the isosurface threshold value "minValue"
77  private float minValue; // threshold value to determine whether a point is in or out of the surface
78  private Vector4[] points;
79 
80  private int depthSize; // dimension of a side
81  private int sliceSize; // dimension of a slice
82  private Vector3 center; // value used to center the objects to be drawn on 0,0,0 on the scene
83 
84  private int nbIndTriangles; // size of mesh.triangles
85  private int nbIndVertices; // size of mesh.vertices
86  private bool[] marchedCubes; // list of length "dataLength", indicates if cubes already have been marched
87  private bool[] marchedCubes2;
88  private bool[] queuedCubes;
89  private bool[] queuedCubes2;
90  private int[] vertIndexGrid; // list of length "dataLength", records indexes of previously created vertices, so we can build the triangle indexes list using unique vertices
91  private int[] vertIndexGrid2;
92  private int dataLength; // nb of points in 2 slices, = ncellsX * ncellsZ * 2
93 // private int fullDataLength; // total nb of points, = ncellsX * ncellsY * ncellsZ
94  private int indFound;
95 
96  private int numVertices; // indice of created vertices
97  private int[] TMPtriangles; // to store mesh.triangles, before we know the size
98  private Vector3[] TMPvertices; // to store mesh.vertices, before we know the size
99  private int[] Mtriangles; // to store final array of mesh.triangles
100  private Vector3[] Mvertices; // to store final array of mesh.vertices
101 
102  private Vector4[] verts;
103  private Vector3[] intVerts;
104  private int cubeIndex; // list cubeIndex of already marched cubes
105  private int[] index;
106  private int edgeIndex;
107  private Vector4 p; // Vector4 returned by LinearInterp, to calculate vertices position
108  private int surfaceNb; // identify each part of the created surface, split by parts of ~64900 vertices
109  private GenerateMesh GMInstance;
110 
111  private bool[] foundVert;
112  private Queue<int> IndQueue;
113  private Queue<int> IndQueue2;
114 // private Queue<int[]> PosQueue;
115 // private Queue<int[]> PosQueue2;
116 
117  private Queue<int> XPosQueue;
118  private Queue<int> YPosQueue;
119  private Queue<int> ZPosQueue;
120  private Queue<int> XPosQueue2;
121  private Queue<int> YPosQueue2;
122  private Queue<int> ZPosQueue2;
123 
124  private int ind;
125  private int fullInd;
126 // private int[] gridPos;
127  private int indice;
128 // private int[] gridTemp;
129  private int sliceNb;
130 
131  private int XPos;
132  private int YPos;
133  private int ZPos;
134 
135  private Vector3 Delta;
136  private Vector3 Origin;
137  private string tag;
138 
139  private Color[] colors;
140 
141  public void MCRecMain(int NX, int NY, int NZ, float minV, Vector4[] P, float tolV,
142  bool reversedThreshold, Vector3 delta, Vector3 origin,
143  Color[] colors = null) {
144 
145 
146  // TODO : Séparer les variables par type d'utilisation (MC, global data, ...)
147 
148  this.colors = colors;
149 
150 // gridTemp = new int[3];
151  // /!\ ORIENTATION OF THE GRID : X, Y classical, Z as depth (may be reversed as : Z,Y,X)
152  ncellsX = NX;
153  ncellsY = NY;
154  ncellsZ = NZ;
155 // toleranceValue = tolV;
156  minValue = minV;
157 
158  foundVert = new bool[12];
159  depthSize = ncellsX; // fix depth axe
160  sliceSize = (ncellsZ) * depthSize;
161 
162  dataLength = ncellsX * ncellsZ * 2;
163 // fullDataLength = ncellsX * ncellsY * ncellsZ;
164 
165  surfaceNb = 0;
166 
167 // center = new Vector3 (-ncellsX/2.0f, -ncellsY/2.0f, -ncellsZ/2.0f); // value used to center the objects to be drawn on 0,0,0 on the scene
168  center = new Vector3(0f,0f,0f);
169 
170  nbIndVertices = 0;
171  nbIndTriangles = 0;
172  numVertices = 0;
173  MarchingModel.TMPvertices = new Vector3[64950]; // unity limit
174  MarchingModel.TMPColors = new Color[64950];
175  MarchingModel.TMPtriangles = new int[65000*6]; // size is arbitrary, it can be changed to less or more depending of the type of data (for triangles)
176  IndQueue = new Queue<int>();
177  IndQueue2 = new Queue<int>();
178 // PosQueue = new Queue<int[]>();
179 // PosQueue2 = new Queue<int[]>();
180 
181  XPosQueue = new Queue<int>();
182  YPosQueue = new Queue<int>();
183  ZPosQueue = new Queue<int>();
184  XPosQueue2 = new Queue<int>();
185  YPosQueue2 = new Queue<int>();
186  ZPosQueue2 = new Queue<int>();
187 
188  verts = new Vector4[8];
189  points = new Vector4[dataLength];
190  intVerts = new Vector3[12];
191  index = new int[3];
192 // gridPos = new int[3];
193  p = new Vector4();
194  vertIndexGrid = new int[sliceSize * 12];
195  vertIndexGrid2 = new int[sliceSize * 12];
196  GMInstance = new GenerateMesh();
197 
198  marchedCubes = new bool[sliceSize]; // initialize to false, no cubes have been marched nor queued
199  marchedCubes2 = new bool[sliceSize];
200  queuedCubes = new bool[sliceSize];
201  queuedCubes2 = new bool[sliceSize];
202 
203 // Debug.Log("MC tables test : "+MarchingCubesTables.edgeTable[1]);
204  Delta = delta;
205  Origin = origin;
206 
207 // Debug.Log("Origin: "+ Origin);
208 
209  sliceNb = ncellsY - 2;
210 
211  for (int i = 0; i < sliceSize; i++) { // TODO : remove useless ones
212  marchedCubes[i] = false;
213  marchedCubes2[i] = false;
214  queuedCubes2[i] = false;
215  queuedCubes[i] = false;
216  }
217  for (int i=0; i < sliceSize * 12; i++) {
218  vertIndexGrid[i] = -1;
219  vertIndexGrid2[i] = -1;
220  }
221  MarchingModel.Mtriangles = new int[nbIndTriangles]; // allocate mesh triangles and vertices, final size
222  MarchingModel.Mvertices = new Vector3[nbIndVertices];
223  MarchingModel.MColors = new Color[nbIndVertices];
224 
225  if (reversedThreshold)
226  MCRecReversedThres(P);
227  else
228  MCRec(P); // begin marching cubes
229 
230  MarchingModel.Mtriangles = new int[nbIndTriangles]; // allocate mesh triangles and vertices, final size
231  MarchingModel.Mvertices = new Vector3[nbIndVertices];
232  MarchingModel.MColors = new Color[nbIndVertices];
233 
234  for (int i=0; i < nbIndTriangles; i++)
235  MarchingModel.Mtriangles[i] = MarchingModel.TMPtriangles[i];
236 
237  for (int i=0; i < nbIndVertices; i++){
238  MarchingModel.Mvertices[i] = MarchingModel.TMPvertices[i];
239  MarchingModel.MColors[i] = MarchingModel.TMPColors[i];
240  }
241  GMInstance.GM(MarchingModel.Mvertices, MarchingModel.Mtriangles, center, surfaceNb,MarchingModel.MColors); // generate display, potentially split on several objects
242 
243  SetAllVariablesToNull() ;
244  }
245 
246  void MCInitAndFindEachIndex (Vector4[] P) { // find all points inside the surface, stored in "indList"
247 
248 // Debug.Log ("MCInitAnd.., sliceNb : "+sliceNb);
249 
250  for (int j = 0; j < dataLength; j++) { // record density values for the current slice
251  fullInd = sliceNb * sliceSize + j;
252  points[j] = P[fullInd];
253  }
254 
255  for (int j = 0; j < sliceSize; j++) {
256  if (points[j].w >= minValue && !queuedCubes[j]) { // queue intersecting cubes
257  if ((int)points[j].x != ncellsX-1 && (int)points[j].y != ncellsY-1 && (int)points[j].z != ncellsZ-1) {
258  IndQueue.Enqueue(j);
259  XPosQueue.Enqueue((int)points[j].x);
260  YPosQueue.Enqueue((int)points[j].y);
261  ZPosQueue.Enqueue((int)points[j].z);
262 // Debug.Log ("MCINit and find queued : "+(int)points[j].x+", "+(int)points[j].y+", "+(int)points[j].z);
263  queuedCubes[j] = true;
264  }
265  }
266  }
267 // Debug.Log ("IndQueue.Count : "+IndQueue.Count);
268  }
269 
270  void MCInitAndFindEachIndex_2 (Vector4[] P) { // find all points inside the surface, stored in "indList"
271 
272 // Debug.Log ("MCInitAnd..2, sliceNb : "+sliceNb);
273 
274  for (int j = 0; j < dataLength; j++) {
275  fullInd = sliceNb * sliceSize + j;
276  points[j] = P[fullInd];
277  }
278 
279  for (int j = 0; j < sliceSize; j++) {
280  if (points[j].w >= minValue && !queuedCubes2[j]) {
281  if ((int)points[j].x != ncellsX-1 && (int)points[j].y != ncellsY-1 && (int)points[j].z != ncellsZ-1) {
282  IndQueue2.Enqueue(j);
283  XPosQueue2.Enqueue((int)points[j].x);
284  YPosQueue2.Enqueue((int)points[j].y);
285  ZPosQueue2.Enqueue((int)points[j].z);
286 // Debug.Log ("MCINit and find 2 queued : "+(int)points[j].x+", "+(int)points[j].y+", "+(int)points[j].z);
287  queuedCubes2[j] = true;
288  }
289  }
290  }
291 // Debug.Log ("IndQueue2.Count : "+IndQueue2.Count);
292  }
293 
294  void MCRec(Vector4[] P) { // proceed recursive marching cubes on yet unproceeded cubes, starting cube at index : "ind"
295 
296  MCInitAndFindEachIndex (P); // find all points inside the surface at first execution
297 
298  while (sliceNb >= 2) {
299 
300 // Debug.Log("IndQueue.Count : "+IndQueue.Count);
301  while (IndQueue.Count != 0) {
302  ind = IndQueue.Dequeue ();
303  XPos = XPosQueue.Dequeue();
304  YPos = YPosQueue.Dequeue();
305  ZPos = ZPosQueue.Dequeue();
306 
307  MarchOneCube();
308 
309  TestFace0 ();
310  TestFace2 ();
311  TestFace3 ();
312  TestFace4 ();
313  TestFace5 ();
314  }
315  sliceNb --;
316 
317  MCInitAndFindEachIndex_2 (P); // find all points inside the surface
318  for (int i = 0; i < sliceSize; i++) {
319  queuedCubes[i] = false;
320  marchedCubes2[i] = false;
321  }
322  for (int i=0; i < sliceSize * 12; i++)
323  vertIndexGrid2[i] = -1;
324 
325 // Debug.Log("IndQueue2.Count : "+IndQueue2.Count);
326  while (IndQueue2.Count != 0) {
327  ind = IndQueue2.Dequeue ();
328  XPos = XPosQueue2.Dequeue();
329  YPos = YPosQueue2.Dequeue();
330  ZPos = ZPosQueue2.Dequeue();
331 
332  MarchOneCube_2();
333 
334  TestFace0_2 ();
335  TestFace2_2 ();
336  TestFace3_2 ();
337  TestFace4_2 ();
338  TestFace5_2 ();
339  }
340  sliceNb --;
341 
342  MCInitAndFindEachIndex (P); // find all points inside the surface
343  for (int i = 0; i < sliceSize; i++) {
344  queuedCubes2[i] = false;
345  marchedCubes[i] = false;
346  }
347  for (int i=0; i < sliceSize * 12; i++)
348  vertIndexGrid[i] = -1;
349  // end of normal cycle
350 
351  if (sliceNb == 0) { // last step
352 
353  while (IndQueue.Count != 0) {
354  ind = IndQueue.Dequeue ();
355  XPos = XPosQueue.Dequeue();
356  YPos = YPosQueue.Dequeue();
357  ZPos = ZPosQueue.Dequeue();
358 
359  MarchOneCube();
360 
361  TestFace0 ();
362  TestFace2 ();
363  TestFace4 ();
364  TestFace5 ();
365  }
366  }
367  }
368  }
369 
370  void TestFace0 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
371  if (XPos > 0 && ((cubeIndex & 1) > 0 || (cubeIndex & 2) > 0 || (cubeIndex & 16) > 0 || (cubeIndex & 32) > 0)) {
372  if (!queuedCubes[ind-1]) {
373 // Debug.Log ("TestFace0 Queued "+((int)points[ind].x-1)+", "+((int)points[ind].y)+", "+(int)points[ind].z);
374  IndQueue.Enqueue(ind - 1);
375  XPosQueue.Enqueue((int)points[ind].x-1);
376  YPosQueue.Enqueue((int)points[ind].y);
377  ZPosQueue.Enqueue((int)points[ind].z);
378  queuedCubes[ind-1] = true;
379  }
380  }
381  }
382  void TestFace0_2 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
383  if (XPos > 0 && ((cubeIndex & 1) > 0 || (cubeIndex & 2) > 0 || (cubeIndex & 16) > 0 || (cubeIndex & 32) > 0)) {
384  if (!queuedCubes2[ind-1]) {
385 // Debug.Log ("TestFace0_2 Queued "+((int)points[ind].x-1)+", "+((int)points[ind].y)+", "+(int)points[ind].z);
386  IndQueue2.Enqueue(ind - 1);
387  XPosQueue2.Enqueue((int)points[ind].x-1);
388  YPosQueue2.Enqueue((int)points[ind].y);
389  ZPosQueue2.Enqueue((int)points[ind].z);
390  queuedCubes2[ind-1] = true;
391  }
392  }
393  }
394  void TestFace2 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
395  if (XPos < ncellsX-2 && ((cubeIndex & 4) > 0 || (cubeIndex & 8) > 0 || (cubeIndex & 64) > 0 || (cubeIndex & 128) > 0)) {
396  if (!queuedCubes[ind+1]) {
397 // Debug.Log ("TestFace2 Queued "+((int)points[ind].x+1)+", "+((int)points[ind].y)+", "+(int)points[ind].z);
398  IndQueue.Enqueue(ind + 1);
399  XPosQueue.Enqueue((int)points[ind].x+1);
400  YPosQueue.Enqueue((int)points[ind].y);
401  ZPosQueue.Enqueue((int)points[ind].z);
402  queuedCubes[ind+1] = true;
403  }
404  }
405  }
406  void TestFace2_2 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
407  if (XPos < ncellsX-2 && ((cubeIndex & 4) > 0 || (cubeIndex & 8) > 0 || (cubeIndex & 64) > 0 || (cubeIndex & 128) > 0)) {
408  if (!queuedCubes2[ind+1]) {
409 // Debug.Log ("TestFace2_2 Queued "+((int)points[ind].x+1)+", "+((int)points[ind].y)+", "+(int)points[ind].z);
410  IndQueue2.Enqueue(ind + 1);
411  XPosQueue2.Enqueue((int)points[ind].x+1);
412  YPosQueue2.Enqueue((int)points[ind].y);
413  ZPosQueue2.Enqueue((int)points[ind].z);
414  queuedCubes2[ind+1] = true;
415  }
416  }
417  }
418  void TestFace3 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
419  // Face 3 : changing slice, invert queues
420  if (YPos > 0 && ((cubeIndex & 1) > 0 || (cubeIndex & 8) > 0 || (cubeIndex & 16) > 0 || (cubeIndex & 128) > 0)) {
421  if (!queuedCubes2[ind]) { // TODO : Useless check ?
422 // Debug.Log ("TestFace3 Queued "+(int)points[ind].x+", "+((int)points[ind].y-1)+", "+((int)points[ind].z));
423  IndQueue2.Enqueue(ind);
424  XPosQueue2.Enqueue((int)points[ind].x);
425  YPosQueue2.Enqueue((int)points[ind].y-1);
426  ZPosQueue2.Enqueue((int)points[ind].z);
427  queuedCubes2[ind] = true;
428  }
429  }
430  }
431  void TestFace3_2 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
432  // Face 3 : changing slice, invert queues
433  if (YPos > 0 && ((cubeIndex & 1) > 0 || (cubeIndex & 8) > 0 || (cubeIndex & 16) > 0 || (cubeIndex & 128) > 0)) {
434  if (!queuedCubes[ind]) { // TODO : Useless check ?
435 // Debug.Log ("TestFace3_2 Queued "+(int)points[ind].x+", "+((int)points[ind].y-1)+", "+((int)points[ind].z));
436  IndQueue.Enqueue(ind);
437  XPosQueue.Enqueue((int)points[ind].x);
438  YPosQueue.Enqueue((int)points[ind].y-1);
439  ZPosQueue.Enqueue((int)points[ind].z);
440  queuedCubes[ind] = true;
441  }
442  }
443  }
444  void TestFace4 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
445  if (ZPos < depthSize-2 && ((cubeIndex & 16) > 0 || (cubeIndex & 32) > 0 || (cubeIndex & 64) > 0 || (cubeIndex & 128) > 0)) {
446  if (!queuedCubes[ind+depthSize]) {
447 // Debug.Log ("TestFace4 Queued "+((int)points[ind].x)+", "+((int)points[ind].y)+", "+((int)points[ind].z+1));
448  IndQueue.Enqueue(ind+depthSize);
449  XPosQueue.Enqueue((int)points[ind].x);
450  YPosQueue.Enqueue((int)points[ind].y);
451  ZPosQueue.Enqueue((int)points[ind].z+1);
452  queuedCubes[ind+depthSize] = true;
453  }
454  }
455  }
456  void TestFace4_2 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
457  if (ZPos < depthSize-2 && ((cubeIndex & 16) > 0 || (cubeIndex & 32) > 0 || (cubeIndex & 64) > 0 || (cubeIndex & 128) > 0)) {
458  if (!queuedCubes2[ind+depthSize]) {
459 // Debug.Log ("TestFace4_2 Queued "+((int)points[ind].x)+", "+((int)points[ind].y-1)+", "+((int)points[ind].z+1));
460  IndQueue2.Enqueue(ind+depthSize);
461  XPosQueue2.Enqueue((int)points[ind].x);
462  YPosQueue2.Enqueue((int)points[ind].y);
463  ZPosQueue2.Enqueue((int)points[ind].z+1);
464  queuedCubes2[ind+depthSize] = true;
465  }
466  }
467  }
468  void TestFace5 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
469  if (ZPos > 0 && ((cubeIndex & 1) > 0 || (cubeIndex & 2) > 0 || (cubeIndex & 4) > 0 || (cubeIndex & 8) > 0)) {
470  if (!queuedCubes[ind-depthSize]) {
471 // Debug.Log ("TestFace5 Queued "+((int)points[ind].x)+", "+((int)points[ind].y)+", "+((int)points[ind].z-1));
472  IndQueue.Enqueue(ind-depthSize);
473  XPosQueue.Enqueue((int)points[ind].x);
474  YPosQueue.Enqueue((int)points[ind].y);
475  ZPosQueue.Enqueue((int)points[ind].z-1);
476  queuedCubes[ind-depthSize] = true;
477  }
478  }
479  }
480  void TestFace5_2 () { // check which adjacent cubes have to be checked depending on previous cube index and grid limits
481  if (ZPos > 0 && ((cubeIndex & 1) > 0 || (cubeIndex & 2) > 0 || (cubeIndex & 4) > 0 || (cubeIndex & 8) > 0)) {
482  if (!queuedCubes2[ind-depthSize]) {
483 // Debug.Log ("TestFace5_2 Queued "+((int)points[ind].x)+", "+((int)points[ind].y)+", "+((int)points[ind].z-1));
484  IndQueue2.Enqueue(ind-depthSize);
485  XPosQueue2.Enqueue((int)points[ind].x);
486  YPosQueue2.Enqueue((int)points[ind].y);
487  ZPosQueue2.Enqueue((int)points[ind].z-1);
488  queuedCubes2[ind-depthSize] = true;
489  }
490  }
491  }
492 
493  void MarchOneCube() { // check how to draw triangles inside current cube, referenced by its indice "ind"
494 
495 // Debug.LogWarning ("MarchOneCube ind : "+ind+", gridPos : "+XPos+", "+YPos+", "+ZPos);
496 
497  if (nbIndVertices > 64900) { // generate several parts of ~64900 vertices if it's going further (Unity limitation at 65000 vertices per GameObject)
498 
499  MarchingModel.Mtriangles = new int[nbIndTriangles]; // allocate mesh triangles and vertices
500  MarchingModel.Mvertices = new Vector3[nbIndVertices];
501  MarchingModel.MColors = new Color[nbIndVertices];
502  for (int i=0; i < nbIndTriangles; i++)
503  MarchingModel.Mtriangles[i] = MarchingModel.TMPtriangles[i];
504 
505  for (int i=0; i < nbIndVertices; i++){
506  MarchingModel.MColors[i] = MarchingModel.TMPColors[i];
507  MarchingModel.Mvertices[i] = MarchingModel.TMPvertices[i];
508  }
509  GMInstance.GM(MarchingModel.Mvertices, MarchingModel.Mtriangles, center, surfaceNb,MarchingModel.MColors); // generate display, eventually split on several objects
510 
511 
512 
513  MarchingModel.TMPtriangles=null;
514  MarchingModel.TMPvertices=null;
515  MarchingModel.TMPColors = null;
516  MarchingModel.TMPtriangles = new int[65000*6]; // reset temporary triangles and vertices
517  MarchingModel.TMPvertices = new Vector3[64950];
518  MarchingModel.TMPColors = new Color[64950];
519  surfaceNb++;
520  nbIndVertices = 0;
521  nbIndTriangles = 0;
522  numVertices = 0;
523  for (int i=0; i < sliceSize * 12; i++){
524  vertIndexGrid[i] = -1;
525  vertIndexGrid2[i] = -1;
526  }
527  }
528 
529  for (int i=0; i < 12; i++)
530  foundVert[i] = false;
531 
532  verts[0] = points[ind]; // list values of current cube's vertex (cube at index : "ind")
533  verts[1] = points[ind + sliceSize];
534  verts[2] = points[ind + sliceSize + 1];
535  verts[3] = points[ind + 1];
536  verts[4] = points[ind + depthSize];
537  verts[5] = points[ind + sliceSize + depthSize];
538  verts[6] = points[ind + sliceSize + depthSize + 1];
539  verts[7] = points[ind + depthSize + 1];
540 
541  cubeIndex = 0;
542  marchedCubes[ind] = true; // make sure this cube is listed as already marched
543 
544 // for (int i=ind*12; i < ind*12+12; i++)
545 // vertIndexGrid[i] = -1;
546 
547  for (int n = 0; n < 8; n++)
548 // if ((verts[n].w >= minValue - toleranceValue) && (verts[n].w <= minValue + toleranceValue))
549  if (verts[n].w >= minValue)
550  cubeIndex |= (1 << n);
551 
552 // Debug.Log("cubeIndex : "+cubeIndex);
553  if (cubeIndex != 0 && cubeIndex != 255) { // if this cube intersects the surface
554 
555  edgeIndex = MarchingCubesTables.edgeTable[cubeIndex];
556 
557  if ((edgeIndex & 1) > 0) {
558  if (ind-1 >= 0) {
559  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+2] != -1) {
560  vertIndexGrid[ind*12] = vertIndexGrid[(ind-1)*12+2];
561  foundVert[0] = true;
562  }
563  }
564  if (ind-depthSize >= 0) {
565  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+4] != -1) {
566  vertIndexGrid[ind*12] = vertIndexGrid[(ind-depthSize)*12+4];
567  foundVert[0] = true;
568  }
569  }
570  if (!foundVert[0]) {
571  LinearInterp(verts[0], verts[1]);
572  intVerts[0] = p;
573  }
574  }
575 
576  if ((edgeIndex & 2) > 0) {
577  if (sliceNb < ncellsY-2) {
578  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+3] != - 1) {
579  vertIndexGrid[ind*12+1] = vertIndexGrid2[ind*12+3];
580  foundVert[1] = true;
581  }
582  }
583  if (ind-depthSize >= 0) {
584  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+5] != -1) {
585  vertIndexGrid[ind*12+1] = vertIndexGrid[(ind-depthSize)*12+5];
586  foundVert[1] = true;
587  }
588  }
589  if (!foundVert[1]) {
590  LinearInterp(verts[1], verts[2]);
591  intVerts[1] = p;
592  }
593  }
594 
595  if ((edgeIndex & 4) > 0) {
596  if (ind+1 <= dataLength) {
597  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12] != -1) {
598  vertIndexGrid[ind*12+2] = vertIndexGrid[(ind+1)*12];
599  foundVert[2] = true;
600  }
601  }
602  if (ind-depthSize >= 0) {
603  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+6] != -1) {
604  vertIndexGrid[ind*12+2] = vertIndexGrid[(ind-depthSize)*12+6];
605  foundVert[2] = true;
606  }
607  }
608  if (!foundVert[2]) {
609  LinearInterp(verts[2], verts[3]);
610  intVerts[2] = p;
611  }
612  }
613 
614  if ((edgeIndex & 8) > 0) {
615  if (ind - depthSize >= 0) {
616  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+7] != -1) {
617  vertIndexGrid[ind*12+3] = vertIndexGrid[(ind-depthSize)*12+7];
618  foundVert[3] = true;
619  }
620  }
621  if (!foundVert[3]) {
622  LinearInterp(verts[3], verts[0]);
623  intVerts[3] = p;
624  }
625  }
626 
627  if ((edgeIndex & 16) > 0) {
628  if (ind-1 >= 0) {
629  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+6] != -1) {
630  vertIndexGrid[ind*12+4] = vertIndexGrid[(ind-1)*12+6];
631  foundVert[4] = true;
632  }
633  }
634  if (ind+depthSize <= dataLength) {
635  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12] != -1) {
636  vertIndexGrid[ind*12+4] = vertIndexGrid[(ind+depthSize)*12];
637  foundVert[4] = true;
638  }
639  }
640  if (!foundVert[4]) {
641  LinearInterp(verts[4], verts[5]);
642  intVerts[4] = p;
643  }
644  }
645 
646  if ((edgeIndex & 32) > 0) {
647  if (sliceNb < ncellsY-2) {
648  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+7] != -1) {
649  vertIndexGrid[ind*12+5] = vertIndexGrid2[ind*12+7];
650  foundVert[5] = true;
651  }
652  }
653  if (ind + depthSize <= dataLength) {
654  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12+1] != -1) {
655  vertIndexGrid[ind*12+5] = vertIndexGrid[(ind+depthSize)*12+1];
656  foundVert[5] = true;
657  }
658  }
659  if (!foundVert[5]) {
660  LinearInterp(verts[5], verts[6]);
661  intVerts[5] = p;
662  }
663  }
664 
665  if ((edgeIndex & 64) > 0) {
666  if (ind+1 <= dataLength) {
667  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12+4] != -1) {
668  vertIndexGrid[ind*12+6] = vertIndexGrid[(ind+1)*12+4];
669  foundVert[6] = true;
670  }
671  }
672  if (ind+ depthSize <= dataLength) {
673  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12+2] != -1) {
674  vertIndexGrid[ind*12+6] = vertIndexGrid[(ind+depthSize)*12+2];
675  foundVert[6] = true;
676  }
677  }
678  if (!foundVert[6]) {
679  LinearInterp(verts[6], verts[7]);
680  intVerts[6] = p;
681  }
682  }
683 
684  if ((edgeIndex & 128) > 0) {
685  if (ind+depthSize <= dataLength) {
686  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12+3] != -1) {
687  vertIndexGrid[ind*12+7] = vertIndexGrid[(ind+depthSize)*12+3];
688  foundVert[7] = true;
689  }
690  }
691  if (!foundVert[7]) {
692  LinearInterp(verts[7], verts[4]);
693  intVerts[7] = p;
694  }
695  }
696 
697  if ((edgeIndex & 256) > 0) {
698  if (ind-1 >= 0) {
699  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+11] != -1) {
700  vertIndexGrid[ind*12+8] = vertIndexGrid[(ind-1)*12+11];
701  foundVert[8] = true;
702  }
703  }
704  if (!foundVert[8]) {
705  LinearInterp(verts[0], verts[4]);
706  intVerts[8] = p;
707  }
708  }
709 
710  if ((edgeIndex & 512) > 0) {
711  if (ind - 1 >= 0) {
712  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+10] != -1) {
713  vertIndexGrid[ind*12+9] = vertIndexGrid[(ind-1)*12+10];
714  foundVert[9] = true;
715  }
716  }
717  if (sliceNb < ncellsY-2) {
718  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+8] != -1) {
719  vertIndexGrid[ind*12+9] = vertIndexGrid2[ind*12+8];
720  foundVert[9] = true;
721  }
722  }
723  if (!foundVert[9]) {
724  LinearInterp(verts[1], verts[5]);
725  intVerts[9] = p;
726  }
727  }
728 
729  if ((edgeIndex & 1024) > 0) {
730  if (ind+1 <= dataLength) {
731  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12+9] != -1) {
732  vertIndexGrid[ind*12+10] = vertIndexGrid[(ind+1)*12+9];
733  foundVert[10] = true;
734  }
735  }
736  if (sliceNb < ncellsY-2) {
737  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+11] != -1) {
738  vertIndexGrid[ind*12+10] = vertIndexGrid2[ind*12+11];
739  foundVert[10] = true;
740  }
741  }
742  if (!foundVert[10]) {
743  LinearInterp(verts[2], verts[6]);
744  intVerts[10] = p;
745  }
746  }
747 
748  if ((edgeIndex & 2048) > 0) {
749  if (ind+1 <= dataLength) {
750  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12+8] != -1) {
751  vertIndexGrid[ind*12+11] = vertIndexGrid[(ind+1)*12+8];
752  foundVert[11] = true;
753  }
754  }
755  if (!foundVert[11]) {
756  LinearInterp(verts[3], verts[7]);
757  intVerts[11] = p;
758  }
759  }
760 
761  for (int n = 0; MarchingCubesTables.triTable[cubeIndex, n] != -1; n += 3) {
762 
763  // TODO : Problem of reversing normals if density threshold < 0
764  index[0] = MarchingCubesTables.triTable[cubeIndex, n];
765  index[1] = MarchingCubesTables.triTable[cubeIndex, n + 1];
766  index[2] = MarchingCubesTables.triTable[cubeIndex, n + 2];
767  for (int h = 0; h < 3; h++)
768  {
769 // Debug.Log("check -1 : "+vertIndexGrid[ind*12+index[h]]);
770  if (vertIndexGrid[ind*12+index[h]] != -1) {
771 // Debug.Log("New triangle index : "+vertIndexGrid[ind*12+index[h]]);
772 // Debug.Log("nbIndTriangles : "+nbIndTriangles);
773  MarchingModel.TMPtriangles[nbIndTriangles] = vertIndexGrid[ind*12+index[h]];
774  }
775  else {
776 // Debug.Log("Triangle index : "+numVertices);
777  MarchingModel.TMPtriangles[nbIndTriangles] = numVertices;
778  vertIndexGrid[ind*12+index[h]] = numVertices;
779  numVertices++;
780  if (Delta.x ==2f){
781  intVerts[index[h]].x -= 18;
782  intVerts[index[h]].y -= 18;
783  intVerts[index[h]].z -= 18;
784  intVerts[index[h]].x /=Delta.x;
785  intVerts[index[h]].y /=Delta.y;
786  intVerts[index[h]].z /=Delta.z;
787  }else{
788  intVerts[index[h]].x *=Delta.x;
789  intVerts[index[h]].y *=Delta.y;
790  intVerts[index[h]].z *=Delta.z;
791  }
792  intVerts[index[h]].x += Origin.x;
793  intVerts[index[h]].y += Origin.y;
794  intVerts[index[h]].z += Origin.z;
795 // Debug.Log("Origin:"+ Origin);
796 
797 
798  // intVerts[index[h]].x += MoleculeModel.Offset.x;//+MoleculeModel.MinValue.x;
799  // intVerts[index[h]].y += MoleculeModel.Offset.y;//+MoleculeModel.MinValue.y;
800  // intVerts[index[h]].z += MoleculeModel.Offset.z;//+MoleculeModel.MinValue.z;
801 // Debug.Log("New vert : "+intVerts[index[h]]);
802  MarchingModel.TMPvertices[nbIndVertices] = intVerts[index[h]];
803 
804  if(colors != null)
805  MarchingModel.TMPColors[nbIndVertices] = colors[sliceNb*sliceSize+ind];
806  else
807  MarchingModel.TMPColors[nbIndVertices] = Color.white;
808  nbIndVertices++;
809  }
810  nbIndTriangles++;
811  }
812  }
813  }
814  }
815 
816  void MarchOneCube_2() { // check how to draw triangles inside current cube, referenced by its indice "ind"
817 
818 // Debug.LogWarning ("MarchOneCube_2 ind : "+ind);
819 
820  if (nbIndVertices > 64900) { // generate several parts of ~64900 vertices if it's going further (Unity limitation at 65000 vertices per GameObject)
821 
822  MarchingModel.Mtriangles = new int[nbIndTriangles]; // allocate mesh triangles and vertices
823  MarchingModel.Mvertices = new Vector3[nbIndVertices];
824  MarchingModel.MColors = new Color[nbIndVertices];
825 
826  for (int i=0; i < nbIndTriangles; i++)
827  MarchingModel.Mtriangles[i] = MarchingModel.TMPtriangles[i];
828 
829  for (int i=0; i < nbIndVertices; i++){
830  MarchingModel.Mvertices[i] = MarchingModel.TMPvertices[i];
831  MarchingModel.MColors[i]=MarchingModel.TMPColors[i];
832  }
833  GMInstance.GM(MarchingModel.Mvertices, MarchingModel.Mtriangles, center, surfaceNb,MarchingModel.MColors); // generate display, eventually split on several objects
834 
835  MarchingModel.TMPtriangles = null; // reset temporary triangles and vertices
836  MarchingModel.TMPvertices = null;
837  MarchingModel.TMPColors = null;
838  MarchingModel.TMPtriangles = new int[65000*6]; // reset temporary triangles and vertices
839  MarchingModel.TMPvertices = new Vector3[64950];
840  MarchingModel.TMPColors = new Color[64950];
841  surfaceNb++;
842  nbIndVertices = 0;
843  nbIndTriangles = 0;
844  numVertices = 0;
845  for (int i=0; i < sliceSize * 12; i++){
846  vertIndexGrid[i] = -1;
847  vertIndexGrid2[i] = -1;
848  }
849  }
850 
851  for (int i=0; i < 12; i++)
852  foundVert[i] = false;
853 
854  verts[0] = points[ind]; // list values of current cube's vertex (cube at index : "ind")
855  verts[1] = points[ind + sliceSize];
856  verts[2] = points[ind + sliceSize + 1];
857  verts[3] = points[ind + 1];
858  verts[4] = points[ind + depthSize];
859  verts[5] = points[ind + sliceSize + depthSize];
860  verts[6] = points[ind + sliceSize + depthSize + 1];
861  verts[7] = points[ind + depthSize + 1];
862 
863  cubeIndex = 0;
864  marchedCubes2[ind] = true; // make sure this cube is listed as already marched
865 
866 // for (int i=ind*12; i < ind*12+12; i++)
867 // vertIndexGrid[i] = -1;
868 
869  for (int n = 0; n < 8; n++)
870 // if ((verts[n].w >= minValue - toleranceValue) && (verts[n].w <= minValue + toleranceValue))
871  if (verts[n].w >= minValue)
872  cubeIndex |= (1 << n);
873 
874 // Debug.Log("cubeIndex : "+cubeIndex);
875  if (cubeIndex != 0 && cubeIndex != 255) { // if this cube intersects the surface
876 
877  edgeIndex = MarchingCubesTables.edgeTable[cubeIndex];
878 
879  if ((edgeIndex & 1) > 0) {
880  if (ind-1 >= 0) {
881  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+2] != -1) {
882  vertIndexGrid2[ind*12] = vertIndexGrid2[(ind-1)*12+2];
883  foundVert[0] = true;
884  }
885  }
886  if (ind-depthSize >= 0) {
887  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+4] != -1) {
888  vertIndexGrid2[ind*12] = vertIndexGrid2[(ind-depthSize)*12+4];
889  foundVert[0] = true;
890  }
891  }
892  if (!foundVert[0]) {
893  LinearInterp(verts[0], verts[1]);
894  intVerts[0] = p;
895  }
896  }
897 
898  if ((edgeIndex & 2) > 0) {
899  if (sliceNb < ncellsY-2) {
900  if (marchedCubes[ind] && vertIndexGrid[ind*12+3] != - 1) {
901  vertIndexGrid2[ind*12+1] = vertIndexGrid[ind*12+3];
902  foundVert[1] = true;
903  }
904  }
905  if (ind-depthSize >= 0) {
906  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+5] != -1) {
907  vertIndexGrid2[ind*12+1] = vertIndexGrid2[(ind-depthSize)*12+5];
908  foundVert[1] = true;
909  }
910  }
911  if (!foundVert[1]) {
912  LinearInterp(verts[1], verts[2]);
913  intVerts[1] = p;
914  }
915  }
916 
917  if ((edgeIndex & 4) > 0) {
918  if (ind+1 <= dataLength) {
919  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12] != -1) {
920  vertIndexGrid2[ind*12+2] = vertIndexGrid2[(ind+1)*12];
921  foundVert[2] = true;
922  }
923  }
924  if (ind-depthSize >= 0) {
925  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+6] != -1) {
926  vertIndexGrid2[ind*12+2] = vertIndexGrid2[(ind-depthSize)*12+6];
927  foundVert[2] = true;
928  }
929  }
930  if (!foundVert[2]) {
931  LinearInterp(verts[2], verts[3]);
932  intVerts[2] = p;
933  }
934  }
935 
936  if ((edgeIndex & 8) > 0) {
937  if (ind - depthSize >= 0) {
938  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+7] != -1) {
939  vertIndexGrid2[ind*12+3] = vertIndexGrid2[(ind-depthSize)*12+7];
940  foundVert[3] = true;
941  }
942  }
943  if (!foundVert[3]) {
944  LinearInterp(verts[3], verts[0]);
945  intVerts[3] = p;
946  }
947  }
948 
949  if ((edgeIndex & 16) > 0) {
950  if (ind-1 >= 0) {
951  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+6] != -1) {
952  vertIndexGrid2[ind*12+4] = vertIndexGrid2[(ind-1)*12+6];
953  foundVert[4] = true;
954  }
955  }
956  if (ind+depthSize <= dataLength) {
957  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12] != -1) {
958  vertIndexGrid2[ind*12+4] = vertIndexGrid2[(ind+depthSize)*12];
959  foundVert[4] = true;
960  }
961  }
962  if (!foundVert[4]) {
963  LinearInterp(verts[4], verts[5]);
964  intVerts[4] = p;
965  }
966  }
967 
968  if ((edgeIndex & 32) > 0) {
969  if (sliceNb < ncellsY-2) {
970  if (marchedCubes[ind] && vertIndexGrid[ind*12+7] != -1) {
971  vertIndexGrid2[ind*12+5] = vertIndexGrid[ind*12+7];
972  foundVert[5] = true;
973  }
974  }
975  if (ind + depthSize <= dataLength) {
976  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12+1] != -1) {
977  vertIndexGrid2[ind*12+5] = vertIndexGrid2[(ind+depthSize)*12+1];
978  foundVert[5] = true;
979  }
980  }
981  if (!foundVert[5]) {
982  LinearInterp(verts[5], verts[6]);
983  intVerts[5] = p;
984  }
985  }
986 
987  if ((edgeIndex & 64) > 0) {
988  if (ind+1 <= dataLength) {
989  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12+4] != -1) {
990  vertIndexGrid2[ind*12+6] = vertIndexGrid2[(ind+1)*12+4];
991  foundVert[6] = true;
992  }
993  }
994  if (ind+ depthSize <= dataLength) {
995  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12+2] != -1) {
996  vertIndexGrid2[ind*12+6] = vertIndexGrid2[(ind+depthSize)*12+2];
997  foundVert[6] = true;
998  }
999  }
1000  if (!foundVert[6]) {
1001  LinearInterp(verts[6], verts[7]);
1002  intVerts[6] = p;
1003  }
1004  }
1005 
1006  if ((edgeIndex & 128) > 0) {
1007  if (ind+depthSize <= dataLength) {
1008  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12+3] != -1) {
1009  vertIndexGrid2[ind*12+7] = vertIndexGrid2[(ind+depthSize)*12+3];
1010  foundVert[7] = true;
1011  }
1012  }
1013  if (!foundVert[7]) {
1014  LinearInterp(verts[7], verts[4]);
1015  intVerts[7] = p;
1016  }
1017  }
1018 
1019  if ((edgeIndex & 256) > 0) {
1020  if (ind-1 >= 0) {
1021  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+11] != -1) {
1022  vertIndexGrid2[ind*12+8] = vertIndexGrid2[(ind-1)*12+11];
1023  foundVert[8] = true;
1024  }
1025  }
1026  if (!foundVert[8]) {
1027  LinearInterp(verts[0], verts[4]);
1028  intVerts[8] = p;
1029  }
1030  }
1031 
1032  if ((edgeIndex & 512) > 0) {
1033  if (ind - 1 >= 0) {
1034  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+10] != -1) {
1035  vertIndexGrid2[ind*12+9] = vertIndexGrid2[(ind-1)*12+10];
1036  foundVert[9] = true;
1037  }
1038  }
1039  if (sliceNb < ncellsY-2) {
1040  if (marchedCubes[ind] && vertIndexGrid[ind*12+8] != -1) {
1041  vertIndexGrid2[ind*12+9] = vertIndexGrid[ind*12+8];
1042  foundVert[9] = true;
1043  }
1044  }
1045  if (!foundVert[9]) {
1046  LinearInterp(verts[1], verts[5]);
1047  intVerts[9] = p;
1048  }
1049  }
1050 
1051  if ((edgeIndex & 1024) > 0) {
1052  if (ind+1 <= dataLength) {
1053  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12+9] != -1) {
1054  vertIndexGrid2[ind*12+10] = vertIndexGrid2[(ind+1)*12+9];
1055  foundVert[10] = true;
1056  }
1057  }
1058  if (sliceNb < ncellsY-2) {
1059  if (marchedCubes[ind] && vertIndexGrid[ind*12+11] != -1) {
1060  vertIndexGrid2[ind*12+10] = vertIndexGrid[ind*12+11];
1061  foundVert[10] = true;
1062  }
1063  }
1064  if (!foundVert[10]) {
1065  LinearInterp(verts[2], verts[6]);
1066  intVerts[10] = p;
1067  }
1068  }
1069 
1070  if ((edgeIndex & 2048) > 0) {
1071  if (ind+1 <= dataLength) {
1072  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12+8] != -1) {
1073  vertIndexGrid2[ind*12+11] = vertIndexGrid2[(ind+1)*12+8];
1074  foundVert[11] = true;
1075  }
1076  }
1077  if (!foundVert[11]) {
1078  LinearInterp(verts[3], verts[7]);
1079  intVerts[11] = p;
1080  }
1081  }
1082 
1083  for (int n = 0; MarchingCubesTables.triTable[cubeIndex, n] != -1; n += 3) {
1084 
1085  // TODO : Problem of reversing normals if density threshold < 0
1086  index[0] = MarchingCubesTables.triTable[cubeIndex, n];
1087  index[1] = MarchingCubesTables.triTable[cubeIndex, n + 1];
1088  index[2] = MarchingCubesTables.triTable[cubeIndex, n + 2];
1089  for (int h = 0; h < 3; h++)
1090  {
1091  if (vertIndexGrid2[ind*12+index[h]] != -1) {
1092 // Debug.Log("New triangle index : "+vertIndexGrid2[ind*12+index[h]]);
1093  MarchingModel.TMPtriangles[nbIndTriangles] = vertIndexGrid2[ind*12+index[h]];
1094  }
1095  else {
1096 // Debug.Log("Triangle index : "+numVertices);
1097  MarchingModel.TMPtriangles[nbIndTriangles] = numVertices;
1098  vertIndexGrid2[ind*12+index[h]] = numVertices;
1099  numVertices++;
1100  if (Delta.x ==2f){
1101  intVerts[index[h]].x -= 18;
1102  intVerts[index[h]].y -= 18;
1103  intVerts[index[h]].z -= 18;
1104  intVerts[index[h]].x /=Delta.x;
1105  intVerts[index[h]].y /=Delta.y;
1106  intVerts[index[h]].z /=Delta.z;
1107  }else{
1108  intVerts[index[h]].x *=Delta.x;
1109  intVerts[index[h]].y *=Delta.y;
1110  intVerts[index[h]].z *=Delta.z;
1111  }
1112  intVerts[index[h]].x += Origin.x;
1113  intVerts[index[h]].y += Origin.y;
1114  intVerts[index[h]].z += Origin.z;
1115 
1116  // intVerts[index[h]].x += MoleculeModel.Offset.x;//+MoleculeModel.MinValue.x;
1117  // intVerts[index[h]].y += MoleculeModel.Offset.y;//+MoleculeModel.MinValue.y;
1118  // intVerts[index[h]].z += MoleculeModel.Offset.z;//+MoleculeModel.MinValue.z;
1119 
1120 // Debug.Log("New vert : "+intVerts[index[h]]);
1121  MarchingModel.TMPvertices[nbIndVertices] = intVerts[index[h]];
1122  if(colors != null)
1123  MarchingModel.TMPColors[nbIndVertices] = colors[sliceNb*sliceSize+ind];
1124  else
1125  MarchingModel.TMPColors[nbIndVertices] = Color.white;
1126 
1127  nbIndVertices++;
1128  }
1129  nbIndTriangles++;
1130  }
1131  }
1132  }
1133  }
1134 
1135  void LinearInterp(Vector4 p1, Vector4 p2) { // linear interpolation between 2 vertex, to locate vertice position
1136 
1137  if (Math.Abs(p1.w - p2.w) > 0.00001) {
1138  p = p1 + (p2 - p1) / (p2.w - p1.w) * (minValue - p1.w);
1139 // p = (p1 + p2) / 2; // no interpolation, each vertex at the middle of each edge, funky results
1140  }
1141  else
1142  p = p1;
1143 // Debug.Log("p : "+p);
1144  }
1145 
1146  void SetAllVariablesToNull () {
1147  points = null;
1148  marchedCubes = null; // list of length "dataLength", indicates if cubes already have been marched
1149  marchedCubes2 = null;
1150  queuedCubes = null;
1151  queuedCubes2 = null;
1152  vertIndexGrid = null; // list of length "dataLength", records indexes of previously created vertices, so we can build the triangle indexes list using unique vertices
1153  vertIndexGrid2 = null;
1154  MarchingModel.TMPtriangles = null; // to store mesh.triangles, before we know the size
1155  MarchingModel.TMPvertices = null; // to store mesh.vertices, before we know the size
1156  MarchingModel.Mtriangles = null; // to store final array of mesh.triangles
1157  MarchingModel.Mvertices = null; // to store final array of mesh.vertices
1158  verts = null;
1159  intVerts = null;
1160  index = null;
1161  GMInstance = null;
1162  foundVert = null;
1163  IndQueue = null;
1164  IndQueue2 = null;
1165 // PosQueue = null;
1166 // PosQueue2 = null;
1167 // gridPos = null;
1168 // gridTemp = null;
1169  XPosQueue = null;
1170  YPosQueue = null;
1171  ZPosQueue = null;
1172  XPosQueue2 = null;
1173  YPosQueue2 = null;
1174  ZPosQueue2 = null;
1175  }
1176 
1177 
1178 
1179 // ------------------------------------------------------------
1180 // Duplicated code to avoid many checks that could slow the surface creation
1181 // Used to inverse threshold for points and cubes selection
1182 // The only differences in the code concern the selection of points
1183 // using minValue
1184 // ------------------------------------------------------------
1185 
1186 
1187 
1188  void MCInitAndFindEachIndex_Rev (Vector4[] P) { // find all points inside the surface, stored in "indList"
1189 
1190 // Debug.Log ("MCInitAnd.., sliceNb : "+sliceNb);
1191 
1192  for (int j = 0; j < dataLength; j++) { // record density values for the current slice
1193  fullInd = sliceNb * sliceSize + j;
1194  points[j] = P[fullInd];
1195  }
1196 
1197  for (int j = 0; j < sliceSize; j++) {
1198  if (points[j].w <= minValue && !queuedCubes[j]) { // queue intersecting cubes
1199  if ((int)points[j].x != ncellsX-1 && (int)points[j].y != ncellsY-1 && (int)points[j].z != ncellsZ-1) {
1200  IndQueue.Enqueue(j);
1201  XPosQueue.Enqueue((int)points[j].x);
1202  YPosQueue.Enqueue((int)points[j].y);
1203  ZPosQueue.Enqueue((int)points[j].z);
1204 // Debug.Log ("MCINit and find queued : "+(int)points[j].x+", "+(int)points[j].y+", "+(int)points[j].z);
1205  queuedCubes[j] = true;
1206  }
1207  }
1208  }
1209 // Debug.Log ("IndQueue.Count : "+IndQueue.Count);
1210  }
1211 
1212  void MCInitAndFindEachIndex_2_Rev (Vector4[] P) { // find all points inside the surface, stored in "indList"
1213 
1214 // Debug.Log ("MCInitAnd..2, sliceNb : "+sliceNb);
1215 
1216  for (int j = 0; j < dataLength; j++) {
1217  fullInd = sliceNb * sliceSize + j;
1218  points[j] = P[fullInd];
1219  }
1220 
1221  for (int j = 0; j < sliceSize; j++) {
1222  if (points[j].w <= minValue && !queuedCubes2[j]) {
1223  if ((int)points[j].x != ncellsX-1 && (int)points[j].y != ncellsY-1 && (int)points[j].z != ncellsZ-1) {
1224  IndQueue2.Enqueue(j);
1225  XPosQueue2.Enqueue((int)points[j].x);
1226  YPosQueue2.Enqueue((int)points[j].y);
1227  ZPosQueue2.Enqueue((int)points[j].z);
1228 // Debug.Log ("MCINit and find 2 queued : "+(int)points[j].x+", "+(int)points[j].y+", "+(int)points[j].z);
1229  queuedCubes2[j] = true;
1230  }
1231  }
1232  }
1233 // Debug.Log ("IndQueue2.Count : "+IndQueue2.Count);
1234  }
1235 
1236  void MCRecReversedThres(Vector4[] P) { // proceed recursive marching cubes on yet unproceeded cubes, starting cube at index : "ind"
1237 
1238  MCInitAndFindEachIndex_Rev (P); // find all points inside the surface at first execution
1239 
1240  while (sliceNb >= 2) {
1241 
1242 // Debug.Log("IndQueue.Count : "+IndQueue.Count);
1243  while (IndQueue.Count != 0) {
1244  ind = IndQueue.Dequeue ();
1245  XPos = XPosQueue.Dequeue();
1246  YPos = YPosQueue.Dequeue();
1247  ZPos = ZPosQueue.Dequeue();
1248 
1249  MarchOneCube_Rev();
1250 
1251  TestFace0 ();
1252  TestFace2 ();
1253  TestFace3 ();
1254  TestFace4 ();
1255  TestFace5 ();
1256  }
1257  sliceNb --;
1258 
1259  MCInitAndFindEachIndex_2_Rev (P); // find all points inside the surface
1260  for (int i = 0; i < sliceSize; i++) {
1261  queuedCubes[i] = false;
1262  marchedCubes2[i] = false;
1263  }
1264  for (int i=0; i < sliceSize * 12; i++)
1265  vertIndexGrid2[i] = -1;
1266 
1267 // Debug.Log("IndQueue2.Count : "+IndQueue2.Count);
1268  while (IndQueue2.Count != 0) {
1269  ind = IndQueue2.Dequeue ();
1270  XPos = XPosQueue2.Dequeue();
1271  YPos = YPosQueue2.Dequeue();
1272  ZPos = ZPosQueue2.Dequeue();
1273 
1274  MarchOneCube_2_Rev();
1275 
1276  TestFace0_2 ();
1277  TestFace2_2 ();
1278  TestFace3_2 ();
1279  TestFace4_2 ();
1280  TestFace5_2 ();
1281  }
1282  sliceNb --;
1283 
1284  MCInitAndFindEachIndex_Rev (P); // find all points inside the surface
1285  for (int i = 0; i < sliceSize; i++) {
1286  queuedCubes2[i] = false;
1287  marchedCubes[i] = false;
1288  }
1289  for (int i=0; i < sliceSize * 12; i++)
1290  vertIndexGrid[i] = -1;
1291  // end of normal cycle
1292 
1293  if (sliceNb == 0) { // last step
1294 
1295  while (IndQueue.Count != 0) {
1296  ind = IndQueue.Dequeue ();
1297  XPos = XPosQueue.Dequeue();
1298  YPos = YPosQueue.Dequeue();
1299  ZPos = ZPosQueue.Dequeue();
1300 
1301  MarchOneCube_Rev();
1302 
1303  TestFace0 ();
1304  TestFace2 ();
1305  TestFace4 ();
1306  TestFace5 ();
1307  }
1308  }
1309  }
1310  }
1311 
1312  void MarchOneCube_Rev() { // check how to draw triangles inside current cube, referenced by its indice "ind"
1313 
1314 // Debug.LogWarning ("MarchOneCube ind : "+ind);
1315 
1316  if (nbIndVertices > 64900) { // generate several parts of ~64900 vertices if it's going further (Unity limitation at 65000 vertices per GameObject)
1317 
1318  MarchingModel.Mtriangles = new int[nbIndTriangles]; // allocate mesh triangles and vertices
1319  MarchingModel.Mvertices = new Vector3[nbIndVertices];
1320  MarchingModel.MColors = new Color[nbIndVertices];
1321  for (int i=0; i < nbIndTriangles; i++)
1322  MarchingModel.Mtriangles[i] = MarchingModel.TMPtriangles[i];
1323 
1324  for (int i=0; i < nbIndVertices; i++){
1325  MarchingModel.MColors[i] = MarchingModel.TMPColors[i];
1326  MarchingModel.Mvertices[i] = MarchingModel.TMPvertices[i];
1327  }
1328  GMInstance.GM(MarchingModel.Mvertices, MarchingModel.Mtriangles, center, surfaceNb,MarchingModel.MColors); // generate display, eventually split on several objects
1329 
1330 
1331 
1332  MarchingModel.TMPtriangles=null;
1333  MarchingModel.TMPvertices=null;
1334  MarchingModel.TMPColors = null;
1335  MarchingModel.TMPtriangles = new int[65000*6]; // reset temporary triangles and vertices
1336  MarchingModel.TMPvertices = new Vector3[64950];
1337  MarchingModel.TMPColors = new Color[64950];
1338  surfaceNb++;
1339  nbIndVertices = 0;
1340  nbIndTriangles = 0;
1341  numVertices = 0;
1342  for (int i=0; i < sliceSize * 12; i++){
1343  vertIndexGrid[i] = -1;
1344  vertIndexGrid2[i] = -1;
1345  }
1346  }
1347 
1348  for (int i=0; i < 12; i++)
1349  foundVert[i] = false;
1350 
1351  verts[0] = points[ind]; // list values of current cube's vertex (cube at index : "ind")
1352  verts[1] = points[ind + sliceSize];
1353  verts[2] = points[ind + sliceSize + 1];
1354  verts[3] = points[ind + 1];
1355  verts[4] = points[ind + depthSize];
1356  verts[5] = points[ind + sliceSize + depthSize];
1357  verts[6] = points[ind + sliceSize + depthSize + 1];
1358  verts[7] = points[ind + depthSize + 1];
1359 
1360  cubeIndex = 0;
1361  marchedCubes[ind] = true; // make sure this cube is listed as already marched
1362 
1363 // for (int i=ind*12; i < ind*12+12; i++)
1364 // vertIndexGrid[i] = -1;
1365 
1366  for (int n = 0; n < 8; n++)
1367 // if ((verts[n].w >= minValue - toleranceValue) && (verts[n].w <= minValue + toleranceValue))
1368  if (verts[n].w <= minValue)
1369  cubeIndex |= (1 << n);
1370 
1371 // Debug.Log("cubeIndex : "+cubeIndex);
1372  if (cubeIndex != 0 && cubeIndex != 255) { // if this cube intersects the surface
1373 
1374  edgeIndex = MarchingCubesTables.edgeTable[cubeIndex];
1375 
1376  if ((edgeIndex & 1) > 0) {
1377  if (ind-1 >= 0) {
1378  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+2] != -1) {
1379  vertIndexGrid[ind*12] = vertIndexGrid[(ind-1)*12+2];
1380  foundVert[0] = true;
1381  }
1382  }
1383  if (ind-depthSize >= 0) {
1384  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+4] != -1) {
1385  vertIndexGrid[ind*12] = vertIndexGrid[(ind-depthSize)*12+4];
1386  foundVert[0] = true;
1387  }
1388  }
1389  if (!foundVert[0]) {
1390  LinearInterp(verts[0], verts[1]);
1391  intVerts[0] = p;
1392  }
1393  }
1394 
1395  if ((edgeIndex & 2) > 0) {
1396  if (sliceNb < ncellsY-2) {
1397  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+3] != - 1) {
1398  vertIndexGrid[ind*12+1] = vertIndexGrid2[ind*12+3];
1399  foundVert[1] = true;
1400  }
1401  }
1402  if (ind-depthSize >= 0) {
1403  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+5] != -1) {
1404  vertIndexGrid[ind*12+1] = vertIndexGrid[(ind-depthSize)*12+5];
1405  foundVert[1] = true;
1406  }
1407  }
1408  if (!foundVert[1]) {
1409  LinearInterp(verts[1], verts[2]);
1410  intVerts[1] = p;
1411  }
1412  }
1413 
1414  if ((edgeIndex & 4) > 0) {
1415  if (ind+1 <= dataLength) {
1416  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12] != -1) {
1417  vertIndexGrid[ind*12+2] = vertIndexGrid[(ind+1)*12];
1418  foundVert[2] = true;
1419  }
1420  }
1421  if (ind-depthSize >= 0) {
1422  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+6] != -1) {
1423  vertIndexGrid[ind*12+2] = vertIndexGrid[(ind-depthSize)*12+6];
1424  foundVert[2] = true;
1425  }
1426  }
1427  if (!foundVert[2]) {
1428  LinearInterp(verts[2], verts[3]);
1429  intVerts[2] = p;
1430  }
1431  }
1432 
1433  if ((edgeIndex & 8) > 0) {
1434  if (ind - depthSize >= 0) {
1435  if (marchedCubes[ind-depthSize] && vertIndexGrid[(ind-depthSize)*12+7] != -1) {
1436  vertIndexGrid[ind*12+3] = vertIndexGrid[(ind-depthSize)*12+7];
1437  foundVert[3] = true;
1438  }
1439  }
1440  if (!foundVert[3]) {
1441  LinearInterp(verts[3], verts[0]);
1442  intVerts[3] = p;
1443  }
1444  }
1445 
1446  if ((edgeIndex & 16) > 0) {
1447  if (ind-1 >= 0) {
1448  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+6] != -1) {
1449  vertIndexGrid[ind*12+4] = vertIndexGrid[(ind-1)*12+6];
1450  foundVert[4] = true;
1451  }
1452  }
1453  if (ind+depthSize <= dataLength) {
1454  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12] != -1) {
1455  vertIndexGrid[ind*12+4] = vertIndexGrid[(ind+depthSize)*12];
1456  foundVert[4] = true;
1457  }
1458  }
1459  if (!foundVert[4]) {
1460  LinearInterp(verts[4], verts[5]);
1461  intVerts[4] = p;
1462  }
1463  }
1464 
1465  if ((edgeIndex & 32) > 0) {
1466  if (sliceNb < ncellsY-2) {
1467  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+7] != -1) {
1468  vertIndexGrid[ind*12+5] = vertIndexGrid2[ind*12+7];
1469  foundVert[5] = true;
1470  }
1471  }
1472  if (ind + depthSize <= dataLength) {
1473  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12+1] != -1) {
1474  vertIndexGrid[ind*12+5] = vertIndexGrid[(ind+depthSize)*12+1];
1475  foundVert[5] = true;
1476  }
1477  }
1478  if (!foundVert[5]) {
1479  LinearInterp(verts[5], verts[6]);
1480  intVerts[5] = p;
1481  }
1482  }
1483 
1484  if ((edgeIndex & 64) > 0) {
1485  if (ind+1 <= dataLength) {
1486  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12+4] != -1) {
1487  vertIndexGrid[ind*12+6] = vertIndexGrid[(ind+1)*12+4];
1488  foundVert[6] = true;
1489  }
1490  }
1491  if (ind+ depthSize <= dataLength) {
1492  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12+2] != -1) {
1493  vertIndexGrid[ind*12+6] = vertIndexGrid[(ind+depthSize)*12+2];
1494  foundVert[6] = true;
1495  }
1496  }
1497  if (!foundVert[6]) {
1498  LinearInterp(verts[6], verts[7]);
1499  intVerts[6] = p;
1500  }
1501  }
1502 
1503  if ((edgeIndex & 128) > 0) {
1504  if (ind+depthSize <= dataLength) {
1505  if (marchedCubes[ind+depthSize] && vertIndexGrid[(ind+depthSize)*12+3] != -1) {
1506  vertIndexGrid[ind*12+7] = vertIndexGrid[(ind+depthSize)*12+3];
1507  foundVert[7] = true;
1508  }
1509  }
1510  if (!foundVert[7]) {
1511  LinearInterp(verts[7], verts[4]);
1512  intVerts[7] = p;
1513  }
1514  }
1515 
1516  if ((edgeIndex & 256) > 0) {
1517  if (ind-1 >= 0) {
1518  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+11] != -1) {
1519  vertIndexGrid[ind*12+8] = vertIndexGrid[(ind-1)*12+11];
1520  foundVert[8] = true;
1521  }
1522  }
1523  if (!foundVert[8]) {
1524  LinearInterp(verts[0], verts[4]);
1525  intVerts[8] = p;
1526  }
1527  }
1528 
1529  if ((edgeIndex & 512) > 0) {
1530  if (ind - 1 >= 0) {
1531  if (marchedCubes[ind-1] && vertIndexGrid[(ind-1)*12+10] != -1) {
1532  vertIndexGrid[ind*12+9] = vertIndexGrid[(ind-1)*12+10];
1533  foundVert[9] = true;
1534  }
1535  }
1536  if (sliceNb < ncellsY-2) {
1537  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+8] != -1) {
1538  vertIndexGrid[ind*12+9] = vertIndexGrid2[ind*12+8];
1539  foundVert[9] = true;
1540  }
1541  }
1542  if (!foundVert[9]) {
1543  LinearInterp(verts[1], verts[5]);
1544  intVerts[9] = p;
1545  }
1546  }
1547 
1548  if ((edgeIndex & 1024) > 0) {
1549  if (ind+1 <= dataLength) {
1550  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12+9] != -1) {
1551  vertIndexGrid[ind*12+10] = vertIndexGrid[(ind+1)*12+9];
1552  foundVert[10] = true;
1553  }
1554  }
1555  if (sliceNb < ncellsY-2) {
1556  if (marchedCubes2[ind] && vertIndexGrid2[ind*12+11] != -1) {
1557  vertIndexGrid[ind*12+10] = vertIndexGrid2[ind*12+11];
1558  foundVert[10] = true;
1559  }
1560  }
1561  if (!foundVert[10]) {
1562  LinearInterp(verts[2], verts[6]);
1563  intVerts[10] = p;
1564  }
1565  }
1566 
1567  if ((edgeIndex & 2048) > 0) {
1568  if (ind+1 <= dataLength) {
1569  if (marchedCubes[ind+1] && vertIndexGrid[(ind+1)*12+8] != -1) {
1570  vertIndexGrid[ind*12+11] = vertIndexGrid[(ind+1)*12+8];
1571  foundVert[11] = true;
1572  }
1573  }
1574  if (!foundVert[11]) {
1575  LinearInterp(verts[3], verts[7]);
1576  intVerts[11] = p;
1577  }
1578  }
1579 
1580  for (int n = 0; MarchingCubesTables.triTable[cubeIndex, n] != -1; n += 3) {
1581 
1582  // TODO : Problem of reversing normals if density threshold < 0
1583  index[0] = MarchingCubesTables.triTable[cubeIndex, n];
1584  index[1] = MarchingCubesTables.triTable[cubeIndex, n + 1];
1585  index[2] = MarchingCubesTables.triTable[cubeIndex, n + 2];
1586  for (int h = 0; h < 3; h++)
1587  {
1588 // Debug.Log("check -1 : "+vertIndexGrid[ind*12+index[h]]);
1589  if (vertIndexGrid[ind*12+index[h]] != -1) {
1590 // Debug.Log("New triangle index : "+vertIndexGrid[ind*12+index[h]]);
1591 // Debug.Log("nbIndTriangles : "+nbIndTriangles);
1592  MarchingModel.TMPtriangles[nbIndTriangles] = vertIndexGrid[ind*12+index[h]];
1593  }
1594  else {
1595 // Debug.Log("Triangle index : "+numVertices);
1596  MarchingModel.TMPtriangles[nbIndTriangles] = numVertices;
1597  vertIndexGrid[ind*12+index[h]] = numVertices;
1598  numVertices++;
1599  if (Delta.x ==2f){
1600  intVerts[index[h]].x -= 18;
1601  intVerts[index[h]].y -= 18;
1602  intVerts[index[h]].z -= 18;
1603  intVerts[index[h]].x /=Delta.x;
1604  intVerts[index[h]].y /=Delta.y;
1605  intVerts[index[h]].z /=Delta.z;
1606  }else{
1607  intVerts[index[h]].x *=Delta.x;
1608  intVerts[index[h]].y *=Delta.y;
1609  intVerts[index[h]].z *=Delta.z;
1610  }
1611  intVerts[index[h]].x += Origin.x;
1612  intVerts[index[h]].y += Origin.y;
1613  intVerts[index[h]].z += Origin.z;
1614 
1615 
1616  // intVerts[index[h]].x += MoleculeModel.Offset.x;//+MoleculeModel.MinValue.x;
1617  // intVerts[index[h]].y += MoleculeModel.Offset.y;//+MoleculeModel.MinValue.y;
1618  // intVerts[index[h]].z += MoleculeModel.Offset.z;//+MoleculeModel.MinValue.z;
1619 // Debug.Log("New vert : "+intVerts[index[h]]);
1620  MarchingModel.TMPvertices[nbIndVertices] = intVerts[index[h]];
1621 
1622  if(colors != null)
1623  MarchingModel.TMPColors[nbIndVertices] = colors[sliceNb*sliceSize+ind];
1624  else
1625  MarchingModel.TMPColors[nbIndVertices] = Color.white;
1626  nbIndVertices++;
1627  }
1628  nbIndTriangles++;
1629  }
1630  }
1631  }
1632  }
1633 
1634  void MarchOneCube_2_Rev() { // check how to draw triangles inside current cube, referenced by its indice "ind"
1635 
1636 // Debug.LogWarning ("MarchOneCube_2 ind : "+ind+", gridPos : "+XPos+", "+YPos+", "+ZPos);
1637 
1638  if (nbIndVertices > 64900) { // generate several parts of ~64900 vertices if it's going further (Unity limitation at 65000 vertices per GameObject)
1639 
1640  MarchingModel.Mtriangles = new int[nbIndTriangles]; // allocate mesh triangles and vertices
1641  MarchingModel.Mvertices = new Vector3[nbIndVertices];
1642  MarchingModel.MColors = new Color[nbIndVertices];
1643  for (int i=0; i < nbIndTriangles; i++)
1644  MarchingModel.Mtriangles[i] = MarchingModel.TMPtriangles[i];
1645 
1646  for (int i=0; i < nbIndVertices; i++){
1647  MarchingModel.Mvertices[i] = MarchingModel.TMPvertices[i];
1648  MarchingModel.MColors[i]=MarchingModel.TMPColors[i];
1649  }
1650  GMInstance.GM(MarchingModel.Mvertices, MarchingModel.Mtriangles, center, surfaceNb,MarchingModel.MColors); // generate display, eventually split on several objects
1651 
1652  MarchingModel.TMPtriangles = null; // reset temporary triangles and vertices
1653  MarchingModel.TMPvertices = null;
1654  MarchingModel.TMPColors = null;
1655  MarchingModel.TMPtriangles = new int[65000*6]; // reset temporary triangles and vertices
1656  MarchingModel.TMPvertices = new Vector3[64950];
1657  MarchingModel.TMPColors = new Color[64950];
1658  surfaceNb++;
1659  nbIndVertices = 0;
1660  nbIndTriangles = 0;
1661  numVertices = 0;
1662  for (int i=0; i < sliceSize * 12; i++) {
1663  vertIndexGrid2[i] = -1;
1664  vertIndexGrid[i] = -1;
1665  }
1666  }
1667 
1668  for (int i=0; i < 12; i++)
1669  foundVert[i] = false;
1670 
1671  verts[0] = points[ind]; // list values of current cube's vertex (cube at index : "ind")
1672  verts[1] = points[ind + sliceSize];
1673  verts[2] = points[ind + sliceSize + 1];
1674  verts[3] = points[ind + 1];
1675  verts[4] = points[ind + depthSize];
1676  verts[5] = points[ind + sliceSize + depthSize];
1677  verts[6] = points[ind + sliceSize + depthSize + 1];
1678  verts[7] = points[ind + depthSize + 1];
1679 
1680  cubeIndex = 0;
1681  marchedCubes2[ind] = true; // make sure this cube is listed as already marched
1682 
1683 // for (int i=ind*12; i < ind*12+12; i++)
1684 // vertIndexGrid[i] = -1;
1685 
1686  for (int n = 0; n < 8; n++)
1687 // if ((verts[n].w >= minValue - toleranceValue) && (verts[n].w <= minValue + toleranceValue))
1688  if (verts[n].w <= minValue)
1689  cubeIndex |= (1 << n);
1690 
1691 // Debug.Log("cubeIndex : "+cubeIndex);
1692  if (cubeIndex != 0 && cubeIndex != 255) { // if this cube intersects the surface
1693 
1694  edgeIndex = MarchingCubesTables.edgeTable[cubeIndex];
1695 
1696  if ((edgeIndex & 1) > 0) {
1697  if (ind-1 >= 0) {
1698  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+2] != -1) {
1699  vertIndexGrid2[ind*12] = vertIndexGrid2[(ind-1)*12+2];
1700  foundVert[0] = true;
1701  }
1702  }
1703  if (ind-depthSize >= 0) {
1704  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+4] != -1) {
1705  vertIndexGrid2[ind*12] = vertIndexGrid2[(ind-depthSize)*12+4];
1706  foundVert[0] = true;
1707  }
1708  }
1709  if (!foundVert[0]) {
1710  LinearInterp(verts[0], verts[1]);
1711  intVerts[0] = p;
1712  }
1713  }
1714 
1715  if ((edgeIndex & 2) > 0) {
1716  if (sliceNb < ncellsY-2) {
1717  if (marchedCubes[ind] && vertIndexGrid[ind*12+3] != - 1) {
1718  vertIndexGrid2[ind*12+1] = vertIndexGrid[ind*12+3];
1719  foundVert[1] = true;
1720  }
1721  }
1722  if (ind-depthSize >= 0) {
1723  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+5] != -1) {
1724  vertIndexGrid2[ind*12+1] = vertIndexGrid2[(ind-depthSize)*12+5];
1725  foundVert[1] = true;
1726  }
1727  }
1728  if (!foundVert[1]) {
1729  LinearInterp(verts[1], verts[2]);
1730  intVerts[1] = p;
1731  }
1732  }
1733 
1734  if ((edgeIndex & 4) > 0) {
1735  if (ind+1 <= dataLength) {
1736  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12] != -1) {
1737  vertIndexGrid2[ind*12+2] = vertIndexGrid2[(ind+1)*12];
1738  foundVert[2] = true;
1739  }
1740  }
1741  if (ind-depthSize >= 0) {
1742  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+6] != -1) {
1743  vertIndexGrid2[ind*12+2] = vertIndexGrid2[(ind-depthSize)*12+6];
1744  foundVert[2] = true;
1745  }
1746  }
1747  if (!foundVert[2]) {
1748  LinearInterp(verts[2], verts[3]);
1749  intVerts[2] = p;
1750  }
1751  }
1752 
1753  if ((edgeIndex & 8) > 0) {
1754  if (ind - depthSize >= 0) {
1755  if (marchedCubes2[ind-depthSize] && vertIndexGrid2[(ind-depthSize)*12+7] != -1) {
1756  vertIndexGrid2[ind*12+3] = vertIndexGrid2[(ind-depthSize)*12+7];
1757  foundVert[3] = true;
1758  }
1759  }
1760  if (!foundVert[3]) {
1761  LinearInterp(verts[3], verts[0]);
1762  intVerts[3] = p;
1763  }
1764  }
1765 
1766  if ((edgeIndex & 16) > 0) {
1767  if (ind-1 >= 0) {
1768  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+6] != -1) {
1769  vertIndexGrid2[ind*12+4] = vertIndexGrid2[(ind-1)*12+6];
1770  foundVert[4] = true;
1771  }
1772  }
1773  if (ind+depthSize <= dataLength) {
1774  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12] != -1) {
1775  vertIndexGrid2[ind*12+4] = vertIndexGrid2[(ind+depthSize)*12];
1776  foundVert[4] = true;
1777  }
1778  }
1779  if (!foundVert[4]) {
1780  LinearInterp(verts[4], verts[5]);
1781  intVerts[4] = p;
1782  }
1783  }
1784 
1785  if ((edgeIndex & 32) > 0) {
1786  if (sliceNb < ncellsY-2) {
1787  if (marchedCubes[ind] && vertIndexGrid[ind*12+7] != -1) {
1788  vertIndexGrid2[ind*12+5] = vertIndexGrid[ind*12+7];
1789  foundVert[5] = true;
1790  }
1791  }
1792  if (ind + depthSize <= dataLength) {
1793  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12+1] != -1) {
1794  vertIndexGrid2[ind*12+5] = vertIndexGrid2[(ind+depthSize)*12+1];
1795  foundVert[5] = true;
1796  }
1797  }
1798  if (!foundVert[5]) {
1799  LinearInterp(verts[5], verts[6]);
1800  intVerts[5] = p;
1801  }
1802  }
1803 
1804  if ((edgeIndex & 64) > 0) {
1805  if (ind+1 <= dataLength) {
1806  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12+4] != -1) {
1807  vertIndexGrid2[ind*12+6] = vertIndexGrid2[(ind+1)*12+4];
1808  foundVert[6] = true;
1809  }
1810  }
1811  if (ind+ depthSize <= dataLength) {
1812  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12+2] != -1) {
1813  vertIndexGrid2[ind*12+6] = vertIndexGrid2[(ind+depthSize)*12+2];
1814  foundVert[6] = true;
1815  }
1816  }
1817  if (!foundVert[6]) {
1818  LinearInterp(verts[6], verts[7]);
1819  intVerts[6] = p;
1820  }
1821  }
1822 
1823  if ((edgeIndex & 128) > 0) {
1824  if (ind+depthSize <= dataLength) {
1825  if (marchedCubes2[ind+depthSize] && vertIndexGrid2[(ind+depthSize)*12+3] != -1) {
1826  vertIndexGrid2[ind*12+7] = vertIndexGrid2[(ind+depthSize)*12+3];
1827  foundVert[7] = true;
1828  }
1829  }
1830  if (!foundVert[7]) {
1831  LinearInterp(verts[7], verts[4]);
1832  intVerts[7] = p;
1833  }
1834  }
1835 
1836  if ((edgeIndex & 256) > 0) {
1837  if (ind-1 >= 0) {
1838  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+11] != -1) {
1839  vertIndexGrid2[ind*12+8] = vertIndexGrid2[(ind-1)*12+11];
1840  foundVert[8] = true;
1841  }
1842  }
1843  if (!foundVert[8]) {
1844  LinearInterp(verts[0], verts[4]);
1845  intVerts[8] = p;
1846  }
1847  }
1848 
1849  if ((edgeIndex & 512) > 0) {
1850  if (ind - 1 >= 0) {
1851  if (marchedCubes2[ind-1] && vertIndexGrid2[(ind-1)*12+10] != -1) {
1852  vertIndexGrid2[ind*12+9] = vertIndexGrid2[(ind-1)*12+10];
1853  foundVert[9] = true;
1854  }
1855  }
1856  if (sliceNb < ncellsY-2) {
1857  if (marchedCubes[ind] && vertIndexGrid[ind*12+8] != -1) {
1858  vertIndexGrid2[ind*12+9] = vertIndexGrid[ind*12+8];
1859  foundVert[9] = true;
1860  }
1861  }
1862  if (!foundVert[9]) {
1863  LinearInterp(verts[1], verts[5]);
1864  intVerts[9] = p;
1865  }
1866  }
1867 
1868  if ((edgeIndex & 1024) > 0) {
1869  if (ind+1 <= dataLength) {
1870  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12+9] != -1) {
1871  vertIndexGrid2[ind*12+10] = vertIndexGrid2[(ind+1)*12+9];
1872  foundVert[10] = true;
1873  }
1874  }
1875  if (sliceNb < ncellsY-2) {
1876  if (marchedCubes[ind] && vertIndexGrid[ind*12+11] != -1) {
1877  vertIndexGrid2[ind*12+10] = vertIndexGrid[ind*12+11];
1878  foundVert[10] = true;
1879  }
1880  }
1881  if (!foundVert[10]) {
1882  LinearInterp(verts[2], verts[6]);
1883  intVerts[10] = p;
1884  }
1885  }
1886 
1887  if ((edgeIndex & 2048) > 0) {
1888  if (ind+1 <= dataLength) {
1889  if (marchedCubes2[ind+1] && vertIndexGrid2[(ind+1)*12+8] != -1) {
1890  vertIndexGrid2[ind*12+11] = vertIndexGrid2[(ind+1)*12+8];
1891  foundVert[11] = true;
1892  }
1893  }
1894  if (!foundVert[11]) {
1895  LinearInterp(verts[3], verts[7]);
1896  intVerts[11] = p;
1897  }
1898  }
1899 
1900  for (int n = 0; MarchingCubesTables.triTable[cubeIndex, n] != -1; n += 3) {
1901 
1902  // TODO : Problem of reversing normals if density threshold < 0
1903  index[0] = MarchingCubesTables.triTable[cubeIndex, n];
1904  index[1] = MarchingCubesTables.triTable[cubeIndex, n + 1];
1905  index[2] = MarchingCubesTables.triTable[cubeIndex, n + 2];
1906  for (int h = 0; h < 3; h++)
1907  {
1908  if (vertIndexGrid2[ind*12+index[h]] != -1) {
1909 // Debug.Log("New triangle index : "+vertIndexGrid2[ind*12+index[h]]);
1910  MarchingModel.TMPtriangles[nbIndTriangles] = vertIndexGrid2[ind*12+index[h]];
1911  }
1912  else {
1913 // Debug.Log("Triangle index : "+numVertices);
1914  MarchingModel.TMPtriangles[nbIndTriangles] = numVertices;
1915  vertIndexGrid2[ind*12+index[h]] = numVertices;
1916  numVertices++;
1917  if (Delta.x ==2f){
1918  intVerts[index[h]].x -= 18;
1919  intVerts[index[h]].y -= 18;
1920  intVerts[index[h]].z -= 18;
1921  intVerts[index[h]].x /=Delta.x;
1922  intVerts[index[h]].y /=Delta.y;
1923  intVerts[index[h]].z /=Delta.z;
1924  }else{
1925  intVerts[index[h]].x *=Delta.x;
1926  intVerts[index[h]].y *=Delta.y;
1927  intVerts[index[h]].z *=Delta.z;
1928  }
1929 
1930  intVerts[index[h]].x += Origin.x;
1931  intVerts[index[h]].y += Origin.y;
1932  intVerts[index[h]].z += Origin.z;
1933 
1934 
1935  // intVerts[index[h]].x += MoleculeModel.Offset.x;//+MoleculeModel.MinValue.x;
1936  // intVerts[index[h]].y += MoleculeModel.Offset.y;//+MoleculeModel.MinValue.y;
1937  // intVerts[index[h]].z += MoleculeModel.Offset.z;//+MoleculeModel.MinValue.z;
1938 
1939 // Debug.Log("New vert : "+intVerts[index[h]]);
1940  MarchingModel.TMPvertices[nbIndVertices] = intVerts[index[h]];
1941  if(colors != null)
1942  MarchingModel.TMPColors[nbIndVertices] = Color.black;//colors[sliceNb*sliceSize+ind];
1943  else
1944  MarchingModel.TMPColors[nbIndVertices] = Color.white;
1945  nbIndVertices++;
1946  }
1947  nbIndTriangles++;
1948  }
1949  }
1950  }
1951  }
1952 
1953 }
MarchingCubesRec.MCInitAndFindEachIndex_2
void MCInitAndFindEachIndex_2(Vector4[] P)
Definition: MarchingCubesRec.cs:270
MarchingCubesRec.XPosQueue2
Queue< int > XPosQueue2
Definition: MarchingCubesRec.cs:120
MarchingCubesRec.MCRec
void MCRec(Vector4[] P)
Definition: MarchingCubesRec.cs:294
GenerateMesh.GM
void GM(Vector3[] Mvertices, int[] Mtriangles, Vector3 center, int surfaceNb, Color[] Colors)
Definition: GenerateMesh.cs:422
MarchingModel.Mvertices
static Vector3[] Mvertices
Definition: MarchingModel.cs:76
MarchingCubesRec.MCRecMain
void MCRecMain(int NX, int NY, int NZ, float minV, Vector4[] P, float tolV, bool reversedThreshold, Vector3 delta, Vector3 origin, Color[] colors=null)
Definition: MarchingCubesRec.cs:141
MarchingModel.TMPvertices
static Vector3[] TMPvertices
Definition: MarchingModel.cs:74
MarchingModel.Mtriangles
static int[] Mtriangles
Definition: MarchingModel.cs:75
MarchingCubesRec.ZPosQueue2
Queue< int > ZPosQueue2
Definition: MarchingCubesRec.cs:122
MarchingCubesRec.MCInitAndFindEachIndex_2_Rev
void MCInitAndFindEachIndex_2_Rev(Vector4[] P)
Definition: MarchingCubesRec.cs:1212
MarchingCubesRec.IndQueue2
Queue< int > IndQueue2
Definition: MarchingCubesRec.cs:113
MarchingCubesRec.MCInitAndFindEachIndex_Rev
void MCInitAndFindEachIndex_Rev(Vector4[] P)
Definition: MarchingCubesRec.cs:1188
MarchingCubesRec.IndQueue
Queue< int > IndQueue
Definition: MarchingCubesRec.cs:112
MarchingCubesRec.MCInitAndFindEachIndex
void MCInitAndFindEachIndex(Vector4[] P)
Definition: MarchingCubesRec.cs:246
MarchingModel.TMPtriangles
static int[] TMPtriangles
Definition: MarchingModel.cs:73
MarchingCubesRec.XPosQueue
Queue< int > XPosQueue
Definition: MarchingCubesRec.cs:117
MarchingModel.MColors
static Color[] MColors
Definition: MarchingModel.cs:78
MarchingModel.TMPColors
static Color[] TMPColors
Definition: MarchingModel.cs:77
MarchingCubesRec.LinearInterp
void LinearInterp(Vector4 p1, Vector4 p2)
Definition: MarchingCubesRec.cs:1135
MarchingCubesRec.GMInstance
GenerateMesh GMInstance
Definition: MarchingCubesRec.cs:109
MarchingCubesRec.ZPosQueue
Queue< int > ZPosQueue
Definition: MarchingCubesRec.cs:119
MarchingCubesRec.YPosQueue
Queue< int > YPosQueue
Definition: MarchingCubesRec.cs:118
MarchingCubesRec.YPosQueue2
Queue< int > YPosQueue2
Definition: MarchingCubesRec.cs:121
MarchingCubesRec.MCRecReversedThres
void MCRecReversedThres(Vector4[] P)
Definition: MarchingCubesRec.cs:1236