present and ask help to inprove a Catmull-Clark algorithm

hello, first for all this is my first post on forums, so please tell me if I mistake the subforum to post it, or some tags, well now, here we go.
About a month ago I had to make procedural rocks, after tried mixing predone "unit" rocks, I ended up trying to create prisms and aplying subsurf to make it more "realistic"
well, there is a very good catmull-clark implementation on the asset store, and free.
but as I saw on the comments, this can't make the subdivision of triangulated meshes, also cant undo the changes done, and, it subdivide the same for all instances.
in short, I had to make a algorithm that fit my particular problem, so I started doing the standard subsurf to quad meshes and a tool to convert prism from tris to quads, but the behaviour wasn't the same than use the catmull clark directly to the tris, (I checked it on blender and the result for the same mesh, but in quads, is less "spiky").
so I finally spanded it to support tris, I always use it on closed meshes (no "open" faces, all faces are drawed) but I comment a modification to calculate correctly on open meshes (not tested).
well, do the tree parts took me a week, and I was also happy to the behaviour, I previously planned to put it there, as free open code, but I ended up busy doing a simple pathfinder (and translate all the comments, names, and descriptions to english was a pain, too xd).

however this is not perfect, in particular two great problems:

1 - scalates very quickly
the algorithm, works, at 4 iterations... acceptable, after 6... took minutes to do it.
if someone knows how to implement quadratic means to speed up it, instead of use a succesion of iterations, this is your case.

2 - leak meshes
well, the only way that I found to do a different mesh for each instance is make a copy, so, it store a mesh when you start to subdivide (line 218), after that manage this mesh as a separate instance, well, the problem is, when you undo, it restart the original mesh, and I have to discard the modified one (line 244).
also I create a mesh instance for the collider, this is because you probably dont need a highly subdivided mesh to the collider (this drop the frames pretty fast) so, you can make a copy to the current instance of subdivision when it fits well to collision(is necessary add a button for that? I found that always used it from code so... seem to be the most usual situation).

so I invite, I you have time, to make the inproves/modifyes that you considere necessary in order to minimize that problems.
also If you see that I'am wrong on something (or something directly don't works xd) please tell me.

Code (CSharp):

1. using System.Collections;
2. using System.Collections.Generic;
3. using UnityEngine;
4. using UnityEditor;
5. public class CatmullClark : MonoBehaviour
6. {
7.     Mesh horMesh;
8.     Mesh CurrMesh;
9.  
10.     MeshFilter mf;
11.  
12.     int steps = 0;
13.  
14.     /// <summary>
15.     /// very specific function, the mesh need to have all faces formed by
16.     /// two coplanar triangles, otherwise the other triangles will be erased.
17.     /// If the mesh have this configuration (like a irregular, or regular quad prism)
18.     /// it will convert all faces to quads, used to testing purposes, and a specific task.
19.     /// The best case to apply subsurf modifier is have a entire premade quad mesh
20.     /// </summary>
21.     public void Cuadriculate(){
22.         if(mf == null){
23.             mf = (MeshFilter) GetComponent<MeshFilter>();
24.             if (mf == null)
25.                 return;
26.         }
27.  
28.         if(CurrMesh == null){
29.             CurrMesh = mf.mesh;
30.         }
31.  
32.         for (int i = 0; i < CurrMesh.subMeshCount; i++)
33.         {
34.             //iterate between different submeshes
35.             Cuadriculate(i);
36.         }
37.     }
38.  
39.     public void Cuadriculate(int meshIndex)
40.     {
41.         if(CurrMesh.GetTopology (meshIndex) != MeshTopology.Triangles)
42.             return;
43.  
44.         Vector3[] currentShape = (CurrMesh).vertices;
45.         int[] triangles = (CurrMesh).GetIndices (meshIndex);
46.         List<int> quads = new List<int>();
47.  
48.         int total = triangles.Length/3;
49.         int[] currCuad;
50.  
51.         for (int i = 0; i < total; i++)
52.         {
53.             for (int j = i+1; j < total; j++)
54.             {
55.                 currCuad = CompareTriangles(triangles, currentShape, i*3, j*3);
56.                 if(currCuad != null){
57.                     for (int h = 0; h < currCuad.Length; h++){
58.                         quads.Add(currCuad[h]);
59.                     }
60.                 }
61.             }
62.         }
63.  
64.         CurrMesh.SetIndices(quads.ToArray(), MeshTopology.Quads, meshIndex);
65.         CurrMesh.RecalculateBounds();
66.         CurrMesh.RecalculateNormals();
67.         CurrMesh.RecalculateTangents();
68.     }
69.  
70.  
71.     /// <summary>
72.     /// Compare two triangles if they have the same normal, and share two vertices
73.     /// </summary>
74.     /// <param name="indexes"> array to points of each triangle (raw of mesh.triangles)</param>
75.     /// <param name="values"> value in 3d space of each point (raw of mesh.vertices)</param>
76.     /// <param name="indexA">starting point on "indexes" to the first triangle</param>
77.     /// <param name="indexB">starting point on "indexes" to the second triangle</param>
78.     /// <returns>a quad generated by the four external vectors, two shared and both limits</returns>
79.     int[] CompareTriangles(int[] indexes, Vector3[] values, int indexA, int indexB){
80.         if (!compareNormal(indexes, values, indexA, indexB))
81.             return null;
82.         int a = -1, b = -1, c = -1;
83.         int d = 0;
84.         if(compare(indexes, indexA, indexB)){
85.             a = indexB;
86.             d = 1;
87.         }else if(compare(indexes, indexA, indexB+1)){
88.             a = indexB+1;
89.             d = 10;
90.         }else if(compare(indexes, indexA, indexB+2)){
91.             a = indexB+2;
92.             d = 100;
93.         }
94.         //----------------------------------------
95.         if(compare(indexes, indexA+1, indexB)){
96.             b = indexB;
97.             d += 1;
98.         }else if(compare(indexes, indexA+1, indexB+1)){
99.             b = indexB+1;
100.             d += 10;
101.         }else if(compare(indexes, indexA+1, indexB+2)){
102.             b = indexB+2;
103.             d += 100;
104.         }
105.         //----------------------------------------
106.         if(compare(indexes, indexA+2, indexB)){
107.             c = indexB;
108.             d += 1;
109.         }else if(compare(indexes, indexA+2, indexB+1)){
110.             c = indexB+1;
111.             d += 10;
112.         }else if(compare(indexes, indexA+2, indexB+2)){
113.             c = indexB+2;
114.             d += 100;
115.         }
116.  
117.         switch (d)
118.         {
119.             case 110:
120.                 d = indexB;
121.                 break;
122.             case 101:
123.                 d = indexB+1;
124.                 break;
125.             case 11:
126.                 d = indexB+2;
127.                 break;
128.             default:
129.                 return null;
130.         }
131.  
132.         //----------------------------------------
133.         int[] ret = new int [4];
134.         if(0 < a){
135.             if (0 < b){
136.                 ret[0] = indexes[a];
137.                 ret[1] = indexes[d];
138.                 ret[2] = indexes[b];
139.                 ret[3] = indexes[indexA+2];
140.             }else if (0 < c){
141.                 ret[0] = indexes[a];
142.                 ret[1] = indexes[indexA+1];
143.                 ret[2] = indexes[c];
144.                 ret[3] = indexes[d];
145.             }else{
146.                 return null;
147.             }
148.         }else if(0 < b){
149.             if(0 < c){
150.                 ret[0] = indexes[indexA];
151.                 ret[1] = indexes[a];
152.                 ret[2] = indexes[b];
153.                 ret[3] = indexes[d];
154.             }else{
155.                 return null;
156.             }
157.         }else{
158.             return null;
159.         }
160.         return ret;
161.     }
162.  
163.     /// <summary>
164.     /// wrapper to compare two array positions
165.     /// </summary>
166.     /// <param name="index"> used to array of triangles</param>
167.     /// <param name="a">array position a</param>
168.     /// <param name="b">array position b</param>
169.     /// <returns></returns>
170.     bool compare(int[] index, int a, int b){
171.         return (index[a] == index[b]);
172.     }
173.  
174.     /// <summary>
175.     /// compare the normals of two triangles, warning, both triangles must be listed on clockwise order
176.     /// </summary>
177.     /// <param name="indexes"> array to points of each triangle (raw of mesh.triangles)</param>
178.     /// <param name="values"> value in 3d space of each point (raw of mesh.vertices)</param>
179.     /// <param name="a">starting point on "indexes" to the first triangle</param>
180.     /// <param name="b">starting point on "indexes" to the second triangle</param>
181.     /// <returns>true if the normals are the same vector, false otherwise</returns>
182.     bool compareNormal(int[] indexes, Vector3[] values, int a, int b){
183.         Vector3 v1 = Normal(values[indexes[a]],  values[indexes[a+1]], values[indexes[a+2]]);
184.         Vector3 v2 = Normal(values[indexes[b]],  values[indexes[b+1]], values[indexes[b+2]]);
185.  
186.         return (Mathf.Approximately(v1.x, v2.x) && Mathf.Approximately(v1.y, v2.y) && Mathf.Approximately(v1.z, v2.z));
187.     }
188.  
189.     /// <summary>
190.     /// calculate the normal of plane generated by tree points, those points must be listed on clockwise order
191.     /// </summary>
192.     /// <param name="a"></param>
193.     /// <param name="b"></param>
194.     /// <param name="c"></param>
195.     /// <returns> the calculated normal</returns>
196.     Vector3 Normal(Vector3 a, Vector3 b, Vector3 c) {
197.         var dir = Vector3.Cross(b - a, c - a);
198.         var norm = Vector3.Normalize(dir);
199.         return norm;
200.     }
201.  
202.     /* ----------------------------------------------------------------------------
203.      * subdivision part
204.      */
205.     /// <summary>
206.     /// main function of subdivision surface
207.     /// </summary>
208.     public void Subdivide(/*int newSteps*/){
209.         if(mf == null){
210.             mf = (MeshFilter) GetComponent<MeshFilter>();
211.             if (mf == null)
212.                 return;
213.         }
214.         if(steps==0 && horMesh == null)//store original mesh at the beggining
215.             horMesh = mf.sharedMesh;
216.  
217.         if(CurrMesh == null){ //store current mesh at the beggining too, from now modifier just change this instance
218.             CurrMesh = mf.mesh;
219.         }
220.  
221.         steps++;
222.         //steps = newSteps;
223.  
224.         for (int i = 0; i < CurrMesh.subMeshCount; i++)
225.         {
226.             //iterate between different submeshes
227.             if(CurrMesh.GetTopology(i) == MeshTopology.Triangles)
228.                 SubdivideTriangles(i);
229.  
230.             if(CurrMesh.GetTopology (i) == MeshTopology.Quads)
231.                 SubdivideQuads(i);
232.         }
233.      
234.     }
235.  
236.     /// <summary>
237.     /// returns to the starting mesh, this function discard previous and allocate memory to the new instance
238.     /// </summary>
239.     public void ClearMesh(){
240.         if (horMesh == null || CurrMesh == null)
241.             return;
242.         CurrMesh.Clear();
243.         mf.sharedMesh = horMesh;
244.         CurrMesh = mf.mesh;
245.         steps = 0;
246.     }
247.  
248.     /// <summary>
249.     /// edit the mesh collider to match to the current subdivision, if there is no mesh collider component, add one
250.     /// note: add a new mesh instance because further modifies on the graphics mesh don't impact on collider
251.     /// </summary>
252.     public void UpdateCollider(){
253.         MeshCollider MColl = GetComponent<MeshCollider> ();
254.         if(MColl){
255.             MColl.sharedMesh = null;
256.             MColl.sharedMesh = mf.mesh;
257.         }else{
258.             MColl = gameObject.AddComponent<MeshCollider> ();
259.             MColl.sharedMesh = mf.mesh;
260.         }
261.     }
262.  
263.     /// <summary>
264.     /// struct used primary by calculate "edge points" contains the original points that generate the edge,
265.     /// the new point and their position in mesh 3d space
266.     /// </summary>
267.     public struct dataVert
268.     {
269.         public int lim1;
270.         public int lim2;
271.         public int vertIndex;
272.         public Vector3 pos;
273.     }
274.     /// <summary>
275.     /// main function of subdivide, on the case that the mesh is formed by quads
276.     /// </summary>
277.     /// <param name="meshIndex">the index of the submesh that need subdivide</param>
278.     public void SubdivideQuads(int meshIndex){
279.         Vector3[] buffvert = (CurrMesh).vertices;
280.         List<int> cuadrados = new List<int>();
281.         List<int> newQuads = new List<int>();
282.         (CurrMesh).GetIndices (cuadrados, meshIndex);
283.         List<Vector3> vertices = new List<Vector3> ();
284.  
285.         int Total = cuadrados.Count/4;
286.         List<dataVert> dict = new List<dataVert> ();
287.         int[] neighbours = new int[buffvert.Length*4];
288.  
289.         for (int i = 0; i < buffvert.Length ; i++)
290.         {
291.             neighbours[i*4] = -1;
292.             neighbours[i*4+1] = -1;
293.             neighbours[i*4+2] = -1;
294.             neighbours[i*4+3] = -1;
295.             vertices.Add(buffvert[i]);
296.         }
297.  
298.         Vector3[] sumVertFaces = new Vector3[buffvert.Length];
299.  
300.         int a,b,c,d,e,f,g,h;
301.         for (int i = 0; i < Total; i++)
302.         {
303.             //add the vertex center of face
304.             a = cuadrados[i*4];
305.             b = cuadrados[i*4+1];
306.             c = cuadrados[i*4+2];
307.             d = cuadrados[i*4+3];
308.             int center = vertices.Count;
309.             vertices.Add((vertices[a] + vertices[b] + vertices[c] + vertices[d])*.25f);
310.             sumVertFaces [a] += vertices[center];
311.             sumVertFaces [b] += vertices[center];
312.             sumVertFaces [c] += vertices[center];
313.             sumVertFaces [d] += vertices[center];
314.             e = addAxisVertex (a, b, vertices[center], dict, vertices);
315.             f = addAxisVertex (b, c, vertices[center], dict, vertices);
316.             g = addAxisVertex (c, d, vertices[center], dict, vertices);
317.             h = addAxisVertex (d, a, vertices[center], dict, vertices);
318.             AddQuad(newQuads,a,e,center,h);
319.             AddQuad(newQuads,e,b,f,center);
320.             AddQuad(newQuads,center, f,c,g);
321.             AddQuad(newQuads, h, center,g,d);
322.         }
323.  
324.      
325.         for (int i = 0; i < buffvert.Length; i++)
326.         {
327.             vertices[i] = AvgNeighboursQuad(i, dict, buffvert, sumVertFaces);
328.         }
329.  
330.         CurrMesh.SetVertices(vertices);
331.         CurrMesh.SetIndices(newQuads.ToArray(), MeshTopology.Quads, meshIndex);
332.         CurrMesh.RecalculateBounds();
333.         CurrMesh.RecalculateNormals();
334.         CurrMesh.RecalculateTangents();
335.     }
336.  
337.     /// <summary>
338.     /// alloc to the struct used in addAxisVertex, just to evade allocate dinamically this
339.     /// </summary>
340.     dataVert temp;
341.  
342.     /// <summary>
343.     /// add, or keep the calculus of a edge point position, this function is complicated.
344.     /// was designed to be called twice:
345.     /// first execution add the new vertex pos to the "dict" array, and start the calculus of the 3d position
346.     /// second execution, when check if this aren't added yet, end the calculus of the vertex 3d position,
347.     /// using the last "vertex face" that still wasn't calculated on the first time.
348.     /// </summary>
349.     /// <param name="a"> original point a of the edge</param>
350.     /// <param name="b"> second point of the edge</param>
351.     /// <param name="center"> vertex face used in this iteration</param>
352.     /// <param name="dict"> array of structs that allocate the vertex</param>
353.     /// <param name="vertices"></param>
354.     /// <returns></returns>
355.     int addAxisVertex (int a, int b, Vector3 center, List<dataVert> dict, List<Vector3> vertices){
356.         for (int i = 0; i < dict.Count; i++)
357.         {
358.             if ((a == dict[i].lim1 && b == dict[i].lim2) ||
359.             (b == dict[i].lim1 && a == dict[i].lim2)){
360.                 //code to the second call
361.                 temp = dict[i];
362.                 // to make it work on open meshes (not all the edges have two neighbour faces use)
363.                 // ret.pos =  vertices[dict[i].vertIndex]*3*.25f;
364.                 vertices[dict[i].vertIndex] += center * .25f;
365.                 temp.pos += center * .25f;
366.                 dict[i] = temp;
367.                 return dict[i].vertIndex;
368.             }
369.         }
370.  
371.         //code to the first call
372.         dataVert ret = new dataVert();
373.         ret.lim1 = a;
374.         ret.lim2 = b;
375.         // to make it work on open meshes (not all the edges have two neighbour faces use)
376.         // ret.pos = (vertices[a]+vertices[b]+center)*.33f;
377.         ret.pos = (vertices[a]+vertices[b]+center)*.25f;
378.         ret.vertIndex = vertices.Count;
379.         vertices.Add(ret.pos);
380.         dict.Add(ret);
381.         return ret.vertIndex;
382.     }
383.  
384.     /// <summary>
385.     /// whapper to add a new quad to the array of quads,
386.     /// to make it work on the mesh you need to allocate this on clockwise order
387.     /// </summary>
388.     /// <param name="cuads">array of quads vertex indexes (raw from mesh.GetIndices if the mesh topology is quads)</param>
389.     /// <param name="a"> point a</param>
390.     /// <param name="b"> point b</param>
391.     /// <param name="c"> point c</param>
392.     /// <param name="d"> point d</param>
393.     void AddQuad(List<int> cuads, int a, int b, int c, int d){
394.         cuads.Add(a);
395.         cuads.Add(b);
396.         cuads.Add(c);
397.         cuads.Add(d);
398.     }
399.  
400.     /// <summary>
401.     /// calculate the new pos on 3d space of the each original vertex point
402.     /// </summary>
403.     /// <param name="vert"> vertex index in the array of vertices </param>
404.     /// <param name="dict"> list of edges + edge points pos </param>
405.     /// <param name="vertices"> positions of the original vertex on 3d space </param>
406.     /// <param name="vertWeights"> sum of all neighbour vertex faces positions, on 3d space  </param>
407.     /// <returns></returns>
408.     Vector3 AvgNeighboursQuad(int vert, List<dataVert> dict, Vector3[] vertices, Vector3[] vertWeights){
409.         Vector3 R = Vector3.zero;
410.         float n = 0f;
411.         for (int i = 0; i < dict.Count; i++)
412.         {
413.             if ((vert == dict[i].lim1) || (vert == dict[i].lim2)){
414.                 R += (vertices[dict[i].lim1]+vertices[dict[i].lim2]) * .5f;
415.                 n += 1f;
416.             }
417.         }
418.      
419.         return (vertWeights[vert]*(1f/n) +2*R*(1f/n) + (n-3)*vertices[vert])*(1f/n);
420.     }
421.  
422.     /// <summary>
423.     /// main function of subdivide, on the case that the mesh is formed by triangles
424.     /// </summary>
425.     /// <param name="meshIndex">the index of the submesh that need subdivide</param>
426.     public void SubdivideTriangles(int meshIndex){
427.         Vector3[] buffvert = (CurrMesh).vertices;
428.         List<int> triangles = new List<int>();
429.         List<int> newQuads = new List<int>();
430.         (CurrMesh).GetIndices (triangles, meshIndex);
431.         List<Vector3> vertices = new List<Vector3> ();
432.  
433.         int Total = triangles.Count/3;
434.         List<dataVert> dict = new List<dataVert> ();
435.         int[] neighbours = new int[buffvert.Length*3];
436.  
437.         for (int i = 0; i < buffvert.Length ; i++)
438.         {
439.             neighbours[i*3] = -1;
440.             neighbours[i*3+1] = -1;
441.             neighbours[i*3+2] = -1;
442.             vertices.Add(buffvert[i]);
443.         }
444.  
445.         Vector3[] sumVertFaces = new Vector3[buffvert.Length];
446.  
447.         int a,b,c,d,e,f;
448.         for (int i = 0; i < Total; i++)
449.         {
450.             //add the vertex center of face
451.             a = triangles[i*3];
452.             b = triangles[i*3+1];
453.             c = triangles[i*3+2];
454.             int center = vertices.Count;
455.             vertices.Add((vertices[a] + vertices[b] + vertices[c])*.33333333f);
456.             sumVertFaces [a] += vertices[center];
457.             sumVertFaces [b] += vertices[center];
458.             sumVertFaces [c] += vertices[center];
459.             d = addAxisVertex (a, b, vertices[center], dict, vertices);
460.             e = addAxisVertex (b, c, vertices[center], dict, vertices);
461.             f = addAxisVertex (c, a, vertices[center], dict, vertices);
462.             AddQuad(newQuads,a,d,center,f);
463.             AddQuad(newQuads,d,b,e,center);
464.             AddQuad(newQuads, f, center, e, c);
465.         }
466.      
467.         for (int i = 0; i < buffvert.Length; i++)
468.         {
469.             vertices[i] = AvgNeighboursTris(i, dict, buffvert, sumVertFaces);
470.         }
471.  
472.         CurrMesh.SetVertices(vertices);
473.         CurrMesh.SetIndices(newQuads.ToArray(), MeshTopology.Quads, meshIndex);
474.         CurrMesh.RecalculateBounds();
475.         CurrMesh.RecalculateNormals();
476.         CurrMesh.RecalculateTangents();
477.     }
478.  
479.     /// <summary>
480.     /// calculate the new pos on 3d space of the each original vertex point
481.     /// </summary>
482.     /// <param name="vert"> vertex index in the array of vertices </param>
483.     /// <param name="dict"> list of edges + edge points pos </param>
484.     /// <param name="vertices"> positions of the original vertex on 3d space </param>
485.     /// <param name="vertWeights"> sum of all neighbour vertex faces positions, on 3d space </param>
486.     /// <returns></returns>
487.     Vector3 AvgNeighboursTris(int vert, List<dataVert> dict, Vector3[] vertices, Vector3[] vertWeights){
488.         Vector3 R = Vector3.zero;
489.         float n = 0f;
490.         for (int i = 0; i < dict.Count; i++)
491.         {
492.             if ((vert == dict[i].lim1) || (vert == dict[i].lim2)){
493.                 R += (vertices[dict[i].lim1] + vertices[dict[i].lim2]) * .5f;
494.                 n += 1f;
495.             }
496.         }
497.      
498.         return (vertWeights[vert]*(1f/n) + 2*R*(1f/n) + (n-3)*vertices[vert])*(1f/n);
499.     }
500.  
501. }
502.  
503.  
504. [CustomEditor(typeof(CatmullClark))]
505. public class CatmullClarkEditor : Editor
506. {
507.     override public void OnInspectorGUI()
508.     {
509.         CatmullClark CCScript = (CatmullClark)target;
510.  
511.         if(GUILayout.Button ("Convert To Quads"))
512.             CCScript.Cuadriculate();
513.         if(GUILayout.Button ("Subdivide"))
514.             CCScript.Subdivide();
515.         if(GUILayout.Button ("Revert"))
516.             CCScript.ClearMesh();
517.     }
518. }