source: GTP/trunk/Lib/Geom/shared/GTGeometry/src/GeoTreeSimplifier.cpp @ 1070

Revision 1070, 46.9 KB checked in by gumbau, 18 years ago (diff)
Line 
1#include "GeoTreeSimplifier.h"
2#include "Leaf.h"
3
4#include <iostream>
5#include <fstream>
6
7using namespace Geometry;
8using namespace std;
9
10const int VISIT_PARENTS_DEEP=1;
11
12class LeafOctree
13{
14public:
15        float left, right, bottom, top, front, back;
16        std::vector<int> leaves; ///< leaves that are too much big to fit in any of the child leaves of this octree leaf
17        LeafOctree* parent;
18        LeafOctree* children[8];
19        bool PointInsideOctree(float x, float y, float z) const {
20                return (x>=left && x<=right && y>=bottom && y<=top && z>=front && z<=back);
21        }
22        bool HasChildren(void) const { return (children[0] && children[1] && children[2] && children[3] && children[4] && children[5] && children[6] && children[7]); }
23};
24
25TreeSimplifier::TreeSimplifier(const Mesh *m, TIPOFUNC upb)
26{
27        objmesh = m;
28        mtreesimpsequence = new Geometry::TreeSimplificationSequence();
29        activeLeaves =  0;
30        countLeaves      =      0;
31
32        // Output mesh
33        mesh    =       new Geometry::Mesh();
34        *mesh   =       *m;
35
36        //      Sets the progress bar.
37        mUPB    =       upb;
38}
39
40TreeSimplifier::~TreeSimplifier()
41{
42        delete  mtreesimpsequence;
43}
44
45// Returns the simplified mesh.
46Geometry::Mesh *TreeSimplifier::GetMesh ()
47{
48        return  mesh;
49}
50
51Geometry::TreeSimplificationSequence *TreeSimplifier::GetSimplificationSequence()
52{
53        return  mtreesimpsequence;
54}
55/*#include <time.h> // necesario para gettickcount
56#include <map>*/
57void TreeSimplifier::Simplify(Real lodfactor,   Index meshLeaves)
58{
59        long    int     totalv;
60        long    int     newleaf;
61        float   diam;
62
63        totalv  =       0;
64
65/*      // DEBUG
66        FILE *f = fopen("timelog.txt","wt");
67        int initime = time(0);
68        fprintf(f,"initime: %d", initime); fflush(f);
69        // DEBUG*/
70       
71        Mesh2Structure(objmesh, meshLeaves);
72       
73        diam = BoundingSphereDiameter();
74
75        // precalculate the octree to speed up the simplification process
76        const int octree_deep = 2;
77        octree=CreateLeafOctree(octree_deep);
78
79        SetCriteriaOptimized(diam);
80
81        int target_face_count = (int)((lodfactor*0.01f) * activeLeaves);
82        int update_each = (activeLeaves - target_face_count)/68;
83//      int prune_each = octree_deep?(activeLeaves - target_face_count)/octree_deep:0;
84       
85//      FILE *fp = fopen("pruning.txt","wt"); // DEBUG
86
87/*      std::map<float,int> hojas_x_crit;
88        for (int i=0; i<countLeaves; i++)
89                hojas_x_crit[Leaves[i].criteria] = i;
90
91        int *leaf_order = new int[countLeaves];
92        std::map<float,int>::iterator kkit = hojas_x_crit.begin();
93        for (int i=0; i<countLeaves; i++, kkit++)
94                leaf_order[i] = kkit->second;*/
95
96//      int erased_leaves = 0;
97        while (activeLeaves > target_face_count)
98        {
99                static int erased_faces_since_last_update = 0;
100                if (mUPB && erased_faces_since_last_update >= update_each)
101                {
102                        erased_faces_since_last_update=0;
103                        mUPB(1);
104                }
105
106/*              static int erased_faces_since_last_prune = 0;
107                if (octree_deep && erased_faces_since_last_prune >= prune_each)
108                {
109                        erased_faces_since_last_prune=0;
110                        int pruneinitime = time(0);
111                        PruneOctree(octree);
112                        int prunefintime = time(0);
113                        // DEBUG                       
114                        fprintf(fp,"pruning at %d activeLeaves [%d (s)]\n", activeLeaves, prunefintime-pruneinitime); fflush(f);
115                        // DEBUG
116                }*/
117
118                erased_faces_since_last_update++;                               
119//              erased_faces_since_last_prune++;
120
121                newleaf = Collapse(diam);
122                SetCriteria2Optimized(diam, newleaf);
123//              erased_leaves++;
124        }
125
126//      fclose(fp);
127
128/*      int fintime = time(0);
129        fprintf(f,"fintime: %d",fintime); fflush(f);
130        fprintf(f,"total: %d (s)",fintime-initime); fflush(f);
131        fclose(f);*/
132
133        BuildOutputMesh(meshLeaves);   
134}
135
136void TreeSimplifier::Mesh2Structure(const Mesh *mesh, Index meshLeaves)
137{
138        size_t          countv=0;
139        long int        pos=0;
140        long int        v1, v2, v3;
141        long int        triangleID=0;
142
143        countv += vertex_count = mesh->mSubMesh[meshLeaves].mVertexBuffer->mVertexCount;       
144        Vertex  = new float[2*countv][3];
145               
146        size_t  update_each = mesh->mSubMesh[meshLeaves].mVertexBuffer->mVertexCount/5;
147       
148        for (unsigned int j=0; j<mesh->mSubMesh[meshLeaves].mVertexBuffer->mVertexCount; j++)
149        {
150                static size_t   ticks_since_last_update = 0;
151                if (mUPB && ticks_since_last_update>=update_each)
152                {
153                        ticks_since_last_update=0;
154                        mUPB(1);
155                }
156                ticks_since_last_update++;
157               
158                Vertex[pos][0] = mesh->mSubMesh[meshLeaves].mVertexBuffer->mPosition[j].x;
159                Vertex[pos][1] = mesh->mSubMesh[meshLeaves].mVertexBuffer->mPosition[j].y;
160                Vertex[pos][2] = mesh->mSubMesh[meshLeaves].mVertexBuffer->mPosition[j].z;
161                pos++;
162        }
163
164        // Calculate the number of leaves
165        countLeaves += (long)mesh->mSubMesh[meshLeaves].mIndexCount;
166        countLeaves=countLeaves/6; // Each leaf is composed of 6 indices
167
168        if (countLeaves > 0)
169                Leaves = new Leaf[2*countLeaves];
170
171        activeLeaves = countLeaves;
172
173        // Insert leaves into the structure
174        pos=0;
175
176        update_each = mesh->mSubMesh[meshLeaves].mIndexCount / 5;
177       
178        // Each leaf is composed of 6 vertices
179        for (unsigned int j=0; j<mesh->mSubMesh[meshLeaves].mIndexCount; j=j+6)
180        {
181                static size_t   ticks_since_last_update = 0;
182                if (mUPB && ticks_since_last_update>=update_each)
183                {
184                        ticks_since_last_update=0;
185                        mUPB(1);
186                }
187                ticks_since_last_update+=6;
188               
189                // first triangle
190                v1=mesh->mSubMesh[meshLeaves].mIndex[j];
191                v2=mesh->mSubMesh[meshLeaves].mIndex[j+1];
192                v3=mesh->mSubMesh[meshLeaves].mIndex[j+2];
193                Leaves[pos].vertsLeaf[0]=v1;
194                Leaves[pos].vertsLeaf[1]=v2;
195                Leaves[pos].vertsLeaf[2]=v3;
196                Leaves[pos].idTriangle[0]=triangleID;
197                Leaves[pos].exists=true;
198                triangleID++;
199
200                // second triangle
201                v3=mesh->mSubMesh[meshLeaves].mIndex[j+5];
202                Leaves[pos].vertsLeaf[3]=v3;
203                Leaves[pos].idTriangle[1]=triangleID;
204                triangleID++;
205
206                CalculateLeafCenter(Leaves[pos]);
207                CalculateLeafNormal(Leaves[pos]);
208                pos++;
209        }
210}
211
212float TreeSimplifier::max(float a,      float b) const
213{
214        if (a>b) return (a);
215        else return(b);
216}
217
218float TreeSimplifier::min(float a,      float b) const
219{
220        if (a>b) return (b);
221        else return(a);
222}
223
224float TreeSimplifier::distan(float x1, float y1, float z1, float x2, float y2, float z2) const
225{
226        return ((x2-x1)*(x2-x1)) + ((y2-y1)*(y2-y1)) + ((z2-z1)*(z2-z1));
227}
228
229// Calculates the center of a leaf
230void TreeSimplifier::CalculateLeafCenter(Leaf &auxleaf)
231{
232        float   max_x;
233        float   max_y;
234        float   max_z;
235        float   min_x;
236        float   min_y;
237        float   min_z;
238       
239        //x1
240        max_x   =       max(max(Vertex[auxleaf.vertsLeaf[0]][0],Vertex[auxleaf.vertsLeaf[1]][0]),
241                                    max(Vertex[auxleaf.vertsLeaf[2]][0],Vertex[auxleaf.vertsLeaf[3]][0]));
242       
243        min_x   =       min(min(Vertex[auxleaf.vertsLeaf[0]][0],Vertex[auxleaf.vertsLeaf[1]][0]),
244                                        min(Vertex[auxleaf.vertsLeaf[2]][0],Vertex[auxleaf.vertsLeaf[3]][0]));
245
246        auxleaf.center[0] = (max_x + min_x)/2;
247
248        //y1
249        max_y   =       max(max(Vertex[auxleaf.vertsLeaf[0]][1],Vertex[auxleaf.vertsLeaf[1]][1]),
250                                        max(Vertex[auxleaf.vertsLeaf[2]][1],Vertex[auxleaf.vertsLeaf[3]][1]));
251       
252        min_y   =       min(min(Vertex[auxleaf.vertsLeaf[0]][1],Vertex[auxleaf.vertsLeaf[1]][1]),
253                                        min(Vertex[auxleaf.vertsLeaf[2]][1],Vertex[auxleaf.vertsLeaf[3]][1]));
254
255        auxleaf.center[1]       =       (max_y + min_y) / 2;
256       
257        //z1
258        max_z   =       max(max(Vertex[auxleaf.vertsLeaf[0]][2],Vertex[auxleaf.vertsLeaf[1]][2]),
259                                        max(Vertex[auxleaf.vertsLeaf[2]][2],Vertex[auxleaf.vertsLeaf[3]][2]));
260       
261        min_z   =       min(min(Vertex[auxleaf.vertsLeaf[0]][2],Vertex[auxleaf.vertsLeaf[1]][2]),
262                                        min(Vertex[auxleaf.vertsLeaf[2]][2],Vertex[auxleaf.vertsLeaf[3]][2]));
263
264        auxleaf.center[2]       =       (max_z + min_z) / 2;
265}
266
267
268//  calculate the normal vector of a leaf
269void TreeSimplifier::CalculateLeafNormal(Leaf &auxleaf)
270{
271        float onex, oney, onez;
272        float twox, twoy, twoz;
273        float threex, threey, threez;
274
275        onex   = Vertex[auxleaf.vertsLeaf[0]][0]; oney   = Vertex[auxleaf.vertsLeaf[0]][1]; onez   = Vertex[auxleaf.vertsLeaf[0]][2];
276        twox   = Vertex[auxleaf.vertsLeaf[1]][0]; twoy   = Vertex[auxleaf.vertsLeaf[1]][1]; twoz   = Vertex[auxleaf.vertsLeaf[1]][2];
277        threex = Vertex[auxleaf.vertsLeaf[2]][0]; threey = Vertex[auxleaf.vertsLeaf[2]][1]; threez = Vertex[auxleaf.vertsLeaf[2]][2];
278       
279        auxleaf.normal[0] = ((twoz-onez)*(threey-oney)) - ((twoy-oney)*(threez-onez));
280        auxleaf.normal[1] = ((twox-onex)*(threez-onez)) - ((threex-onex)*(twoz-onez));
281        auxleaf.normal[2] = ((threex-onex)*(twoy-oney)) - ((twox-onex)*(threey-oney)); 
282}
283
284
285//  calculate the Hausdorff distance (distance between point clouds)
286/*float TreeSimplifier::Hausdorff(Hoja &leaf1, Hoja& leaf2) const
287{
288        float onex, oney, onez;
289        float twox, twoy, twoz;
290        float threex, threey, threez;
291        float fourx, foury, fourz;
292        float x1, y1, z1;
293        float x2, y2, z2;
294        float x3, y3, z3;
295        float x4, y4, z4;
296        float dist1, dist2, dist3, dist4, distmp, dista, distb, dist;
297
298        onex=Vertex[leaf1.vertsLeaf[0]][0]; oney= Vertex[leaf1.vertsLeaf[0]][1]; onez = Vertex[leaf1.vertsLeaf[0]][2];
299        twox = Vertex[leaf1.vertsLeaf[1]][0]; twoy = Vertex[leaf1.vertsLeaf[1]][1]; twoz = Vertex[leaf1.vertsLeaf[1]][2];
300        threex = Vertex[leaf1.vertsLeaf[2]][0]; threey = Vertex[leaf1.vertsLeaf[2]][1]; threez = Vertex[leaf1.vertsLeaf[2]][2];
301        fourx = Vertex[leaf1.vertsLeaf[3]][0]; foury = Vertex[leaf1.vertsLeaf[3]][1]; fourz = Vertex[leaf1.vertsLeaf[3]][2];
302       
303
304        x1 = Vertex[leaf2.vertsLeaf[0]][0]; y1 = Vertex[leaf2.vertsLeaf[0]][1]; z1 = Vertex[leaf2.vertsLeaf[0]][2];
305        x2 = Vertex[leaf2.vertsLeaf[1]][0]; y2 = Vertex[leaf2.vertsLeaf[1]][1]; z2 = Vertex[leaf2.vertsLeaf[1]][2];
306        x3 = Vertex[leaf2.vertsLeaf[2]][0]; y3 = Vertex[leaf2.vertsLeaf[2]][1]; z3 = Vertex[leaf2.vertsLeaf[2]][2];
307        x4 = Vertex[leaf2.vertsLeaf[3]][0]; y4 = Vertex[leaf2.vertsLeaf[3]][1]; z4 = Vertex[leaf2.vertsLeaf[3]][2];
308
309        dist1 = distan ( onex, oney,onez,x1,y1,z1);
310
311        distmp = distan ( onex, oney,onez,x2,y2,z2);
312        if ( distmp < dist1) dist1 = distmp;
313
314        distmp = distan ( onex, oney,onez,x3,y3,z3);
315        if ( distmp < dist1) dist1 = distmp;
316
317        distmp = distan ( onex, oney,onez,x4,y4,z4);
318        if ( distmp < dist1) dist1 = distmp;
319
320        dist2 = distan ( twox, twoy,twoz,x1,y1,z1);
321       
322        distmp = distan ( twox, twoy,twoz,x2,y2,z2);
323        if ( distmp < dist2) dist2 = distmp;
324       
325        distmp = distan ( twox, twoy,twoz,x3,y3,z3);
326        if ( distmp < dist2) dist2 = distmp;
327       
328        distmp = distan ( twox, twoy,twoz,x4,y4,z4);
329        if ( distmp < dist2) dist2 = distmp;
330
331
332        dist3 = distan ( threex, threey,threez,x1,y1,z1);
333       
334        distmp = distan ( threex, threey,threez,x2,y2,z2);
335        if ( distmp < dist3) dist3 = distmp;
336
337        distmp = distan ( threex, threey,threez,x3,y3,z3);
338        if ( distmp < dist3) dist3 = distmp;
339
340        distmp = distan ( threex, threey,threez,x4,y4,z4);
341        if ( distmp < dist3) dist3 = distmp;
342
343
344
345        dist4 = distan ( fourx, foury,fourz,x1,y1,z1);
346       
347        distmp = distan ( fourx, foury,fourz,x2,y2,z2);
348        if ( distmp < dist4) dist4 = distmp;
349
350        distmp = distan ( fourx, foury,fourz,x3,y3,z3);
351        if ( distmp < dist4) dist4 = distmp;
352       
353        distmp = distan ( fourx, foury,fourz,x4,y4,z4);
354        if ( distmp < dist4) dist4 = distmp;
355
356
357        dista =  max(dist1, max(dist2, max(dist3, dist4)));
358
359        dist1 = distan ( x1,y1,z1, onex, oney, onez);
360
361        distmp = distan ( x1,y1,z1, twox, twoy, twoz);
362        if ( distmp < dist1) dist1 = distmp;
363
364        distmp = distan ( x1,y1,z1, threex, threey, threez);
365        if ( distmp < dist1) dist1 = distmp;
366
367        distmp = distan ( x1,y1,z1, fourx, foury, fourz);
368        if ( distmp < dist1) dist1 = distmp;
369
370        dist2 = distan ( x2,y2,z2, onex, oney, onez);
371
372        distmp = distan ( x2,y2,z2, twox, twoy, twoz);
373        if ( distmp < dist2) dist2 = distmp;
374
375        distmp = distan ( x2,y2,z2, threex, threey, threez);
376        if ( distmp < dist2) dist2 = distmp;
377
378        distmp = distan ( x2,y2,z2, fourx, foury, fourz);
379        if ( distmp < dist2) dist2 = distmp;
380
381
382                //3
383        dist3 = distan ( x3,y3,z3, onex, oney, onez);
384
385        distmp = distan ( x3,y3,z3, twox, twoy, twoz);
386        if ( distmp < dist3) dist3 = distmp;
387
388        distmp = distan ( x3,y3,z3, threex, threey, threez);
389        if ( distmp < dist3) dist3 = distmp;
390
391        distmp = distan ( x3,y3,z3, fourx, foury, fourz);
392        if ( distmp < dist3) dist3 = distmp;
393
394                //4
395        dist4 = distan ( x4,y4,z4, onex, oney, onez);
396
397        distmp = distan ( x4,y4,z4, twox, twoy, twoz);
398        if ( distmp < dist4) dist4 = distmp;
399
400        distmp = distan ( x4,y4,z4, threex, threey, threez);
401        if ( distmp < dist4) dist4 = distmp;
402
403        distmp = distan ( x4,y4,z4, fourx, foury, fourz);
404        if ( distmp < dist4) dist4 = distmp;
405
406        //
407        distb =  max(dist1, max(dist2, max(dist3, dist4)));
408
409        dist  = max ( dista, distb);
410        return ( dist);
411       
412}
413*/
414
415//  Calculate the Hausdorff distance (distance between point clouds) (Optimized)
416float TreeSimplifier::HausdorffOptimized(const Leaf &leaf1, const Leaf& leaf2) const
417{
418        float onex, oney, onez;
419        float twox, twoy, twoz;
420        float threex, threey, threez;
421        float fourx, foury, fourz;
422        float x1, y1, z1;
423        float x2, y2, z2;
424        float x3, y3, z3;
425        float x4, y4, z4;
426        float dist1, dist2, dist3, dist4, distmp, dista, distb, dist;
427
428        onex   = Vertex[leaf1.vertsLeaf[0]][0]; oney   = Vertex[leaf1.vertsLeaf[0]][1]; onez   = Vertex[leaf1.vertsLeaf[0]][2];
429        twox   = Vertex[leaf1.vertsLeaf[1]][0]; twoy   = Vertex[leaf1.vertsLeaf[1]][1]; twoz   = Vertex[leaf1.vertsLeaf[1]][2];
430        threex = Vertex[leaf1.vertsLeaf[2]][0]; threey = Vertex[leaf1.vertsLeaf[2]][1]; threez = Vertex[leaf1.vertsLeaf[2]][2];
431        fourx  = Vertex[leaf1.vertsLeaf[3]][0]; foury  = Vertex[leaf1.vertsLeaf[3]][1]; fourz  = Vertex[leaf1.vertsLeaf[3]][2];
432
433        x1 = Vertex[leaf2.vertsLeaf[0]][0]; y1 = Vertex[leaf2.vertsLeaf[0]][1]; z1 = Vertex[leaf2.vertsLeaf[0]][2];
434        x2 = Vertex[leaf2.vertsLeaf[1]][0]; y2 = Vertex[leaf2.vertsLeaf[1]][1]; z2 = Vertex[leaf2.vertsLeaf[1]][2];
435        x3 = Vertex[leaf2.vertsLeaf[2]][0]; y3 = Vertex[leaf2.vertsLeaf[2]][1]; z3 = Vertex[leaf2.vertsLeaf[2]][2];
436        x4 = Vertex[leaf2.vertsLeaf[3]][0]; y4 = Vertex[leaf2.vertsLeaf[3]][1]; z4 = Vertex[leaf2.vertsLeaf[3]][2];
437
438        // variables used to cache distances
439        float one1,one2,one3,one4, two1,two2,two3,two4, three1,three2,three3,three4, four1,four2,four3,four4;
440
441        // Store the minimal distances from each of the 4 vertices to the other 4 vertices
442        one1 = dist1 = distan ( onex, oney,onez,x1,y1,z1);
443       
444        one2 = distmp = distan ( onex, oney,onez,x2,y2,z2);
445        if ( distmp < dist1) dist1 = distmp;
446
447        one3 = distmp = distan ( onex, oney,onez,x3,y3,z3);
448        if ( distmp < dist1) dist1 = distmp;
449
450        one4 = distmp = distan ( onex, oney,onez,x4,y4,z4);
451        if ( distmp < dist1) dist1 = distmp;
452
453       
454        two1 = dist2 = distan ( twox, twoy,twoz,x1,y1,z1);
455       
456        two2 = distmp = distan ( twox, twoy,twoz,x2,y2,z2);
457        if ( distmp < dist2) dist2 = distmp;
458       
459        two3 = distmp = distan ( twox, twoy,twoz,x3,y3,z3);
460        if ( distmp < dist2) dist2 = distmp;
461       
462        two4 = distmp = distan ( twox, twoy,twoz,x4,y4,z4);
463        if ( distmp < dist2) dist2 = distmp;
464
465       
466        three1 = dist3 = distan ( threex, threey,threez,x1,y1,z1);
467       
468        three2 = distmp = distan ( threex, threey,threez,x2,y2,z2);
469        if ( distmp < dist3) dist3 = distmp;
470
471        three3 = distmp = distan ( threex, threey,threez,x3,y3,z3);
472        if ( distmp < dist3) dist3 = distmp;
473
474        three4 = distmp = distan ( threex, threey,threez,x4,y4,z4);
475        if ( distmp < dist3) dist3 = distmp;
476
477
478        four1 = dist4 = distan ( fourx, foury,fourz,x1,y1,z1);
479       
480        four2 = distmp = distan ( fourx, foury,fourz,x2,y2,z2);
481        if ( distmp < dist4) dist4 = distmp;
482
483        four3 = distmp = distan ( fourx, foury,fourz,x3,y3,z3);
484        if ( distmp < dist4) dist4 = distmp;
485       
486        four4 = distmp = distan ( fourx, foury,fourz,x4,y4,z4);
487        if ( distmp < dist4) dist4 = distmp;
488
489        // pick up the maximum value from those 4
490
491        dista =  max(dist1, max(dist2, max(dist3, dist4)));
492
493        // now use the cached distances to perform the reverse distances
494
495        dist1 = one1;  //distan ( x1,y1,z1, onex, oney, onez);
496
497        distmp = two1; //distan ( x1,y1,z1, twox, twoy, twoz);
498        if ( distmp < dist1) dist1 = distmp;
499
500        distmp =  three1; //distan ( x1,y1,z1, threex, threey, threez);
501        if ( distmp < dist1) dist1 = distmp;
502
503        distmp = four1; //distan ( x1,y1,z1, fourx, foury, fourz);
504        if ( distmp < dist1) dist1 = distmp;
505
506        //2
507        dist2 = one2; //distan ( x2,y2,z2, onex, oney, onez);
508
509        distmp = two2; //distan ( x2,y2,z2, twox, twoy, twoz);
510        if ( distmp < dist2) dist2 = distmp;
511
512        distmp = three2; //distan ( x2,y2,z2, threex, threey, threez);
513        if ( distmp < dist2) dist2 = distmp;
514
515        distmp = four2; //distan ( x2,y2,z2, fourx, foury, fourz);
516        if ( distmp < dist2) dist2 = distmp;
517
518
519                //3
520        dist3 = one3; //distan ( x3,y3,z3, onex, oney, onez);
521
522        distmp = two3; //distan ( x3,y3,z3, twox, twoy, twoz);
523        if ( distmp < dist3) dist3 = distmp;
524
525        distmp = three3; //distan ( x3,y3,z3, threex, threey, threez);
526        if ( distmp < dist3) dist3 = distmp;
527
528        distmp = four3; //distan ( x3,y3,z3, fourx, foury, fourz);
529        if ( distmp < dist3) dist3 = distmp;
530
531                //4
532        dist4 = one4; //distan ( x4,y4,z4, onex, oney, onez);
533
534        distmp = two4; //distan ( x4,y4,z4, twox, twoy, twoz);
535        if ( distmp < dist4) dist4 = distmp;
536
537        distmp = three4; //distan ( x4,y4,z4, threex, threey, threez);
538        if ( distmp < dist4) dist4 = distmp;
539
540        distmp = four4; //distan ( x4,y4,z4, fourx, foury, fourz);
541        if ( distmp < dist4) dist4 = distmp;
542
543        //
544        distb =  max(dist1, max(dist2, max(dist3, dist4)));
545
546        dist  = max ( dista, distb);
547        return ( dist);
548       
549}
550
551
552// Calculate the distance between leaves
553/*
554void TreeSimplifier::DistanciaEntreHojas(void)
555{
556        float dist, distmp ;
557        int j;
558
559        for ( int i=0;i<countLeaves;i++)
560        {
561                if ( Leaves[i].existe == true) // si la hoja aun existe
562                {
563                        // inicializo para poder comparar
564                        if ( i == (countLeaves-1)) j = 0;
565                                else j = i+1;
566                        while (Leaves[j].existe == false) j++;
567                        //dist = Leaves[i].Distancia(Leaves[j]);
568                        dist = HausdorffOptimizado( Leaves[i], Leaves[j]);
569
570                        Leaves[i].hoja_cerca = j;
571                        // empiezo los calculos
572                        for ( j =0; j<(countLeaves-1);j++)
573                        {
574                                if ( j == i)
575                                        break;
576                                if ( Leaves[j].existe == true) // si la hoja aun existe
577                                {
578                                        distmp = HausdorffOptimizado(Leaves[i], Leaves[j]);
579                                        if ( distmp < dist )
580                                        {
581                                                dist = distmp;
582                                                Leaves[i].hoja_cerca = j;
583                                        }
584                                }
585                               
586                        }
587                        Leaves[i].dist = dist;
588                }
589        }
590
591}*/
592
593// Returns the leaf which distance is the lowest
594long int TreeSimplifier::MinDistance(void)
595{
596        float           mindist=0.0f;
597        long int        which=-1;
598        int                     i=0;
599
600/*      //init
601        while (Leaves[i].existe != true)
602                i++;
603       
604        mindist = Leaves[i].dist;
605        which   = i;
606*/     
607        // search the minimum
608        for (i = 0; i < countLeaves; i++)
609        {
610                if (Leaves[i].exists)
611                {
612                        if (mindist > Leaves[i].dist || which==-1)
613                        {
614                                mindist = Leaves[i].dist;
615                                which = i;
616                        }
617                }
618
619        }
620        return which;
621}
622
623// Calculate the coplanarity between leaves
624void TreeSimplifier::CoplanarBetweenLeaves(void)
625{
626        float   cop;
627        float   coptmp;
628        int     i;
629        int     j;
630
631        for (i=0;i<countLeaves;i++)
632        {
633                if (Leaves[i].exists)
634                {
635                        if (i == countLeaves-1)
636                                j       =       0;
637                        else
638                                j =     i + 1;
639                       
640                        while (Leaves[j].exists == false)
641                                j++;
642                       
643                        cop     = Leaves[i].Coplanarity(Leaves[j]);
644                        Leaves[i].leafCop       =       j;
645                       
646                        for (j=0; j<countLeaves-1; j++)
647                        {
648                                if (( j != i) && (Leaves[j].exists))
649                                {
650                                        coptmp  =       Leaves[i].Coplanarity(Leaves[j]);
651
652                                        // Take the most coplanar: close to 1
653                                        if (coptmp > cop)
654                                        {
655                                                cop     = coptmp;
656                                                Leaves[i].leafCop = j;
657                                        }
658                                }
659                        }
660                        Leaves[i].coplanar = 1 - cop;
661                }
662        }
663}
664
665void TreeSimplifier::TwoGreater(float *maj, int *indices)
666{
667        float   m1;
668        int             i;
669
670        if      (maj[0] < maj[1])
671        {
672                m1      = maj[0];
673                maj[0] = maj[1];
674                maj[1] = m1;
675
676                i =     indices[0];
677                indices[0] = indices[1];
678                indices[1] = i;
679        }
680
681        if (maj[2] < maj[3])
682        {
683                m1      = maj[2];
684                maj[2] = maj[3];
685                maj[3]  = m1;
686                i       = indices[2];
687                indices[2] = indices[3];
688                indices[3] = i;
689        }
690
691        if  (maj[0] < maj[2])
692        {
693                m1              =       maj[0];
694                maj[0]  =       maj[2];
695                maj[2]  =       m1;
696                i               =       indices[0];
697                indices[0]      =       indices[2];
698                indices[2]      =       i;
699        }
700
701        if (maj[2] < maj[3])
702        {
703                m1              =       maj[2];
704                maj[2]  =       maj[3];
705                maj [3] =       m1;
706
707                i                       =       indices[2];
708                indices[2]      =       indices[3];
709                indices[3]      =       i;
710        }
711
712        if (maj [1] < maj [2])
713        {
714                m1              =       maj[1];
715                maj[1]  =       maj[2];
716                maj [2] =       m1;
717
718                i                       = indices[1];
719                indices[1]      = indices[2];
720                indices[2]      = i;
721        }               
722}
723
724float TreeSimplifier::BoundingSphereDiameter(void) const
725{
726        float xmax, xmin, ymax, ymin, zmax, zmin;
727        float diameter, cx, cy, cz;
728        int j;
729
730        // initialize all vertices to the first
731        xmax = Vertex[Leaves[0].vertsLeaf[0]][0];
732        xmin = xmax;
733
734        ymax = Vertex[Leaves[0].vertsLeaf[0]][1];
735        ymin = ymax;
736
737        zmax = Vertex[Leaves[0].vertsLeaf[0]][2];
738        zmin = zmax;
739
740        // search for max and mins
741
742        int update_each = countLeaves / 5;
743        for (int i = 1; i < countLeaves; i++)
744        {
745                static int ticks_since_last_update = 0;
746                if (mUPB && ticks_since_last_update>=update_each)
747                {
748                        ticks_since_last_update=0;
749                        mUPB(1);
750                }
751                ticks_since_last_update++;
752               
753                for (j=0;j<4;j++)
754                {                       
755                        if ( xmax < Vertex[Leaves[i].vertsLeaf[j]][0]) xmax = Vertex[Leaves[i].vertsLeaf[j]][0];
756                        if ( xmin > Vertex[Leaves[i].vertsLeaf[j]][0]) xmin = Vertex[Leaves[i].vertsLeaf[j]][0];
757
758                        if ( ymax < Vertex[Leaves[i].vertsLeaf[j]][1]) ymax = Vertex[Leaves[i].vertsLeaf[j]][1];
759                        if ( ymin > Vertex[Leaves[i].vertsLeaf[j]][1]) ymin = Vertex[Leaves[i].vertsLeaf[j]][1];
760                       
761                        if ( zmax < Vertex[Leaves[i].vertsLeaf[j]][2]) zmax = Vertex[Leaves[i].vertsLeaf[j]][2];
762                        if ( zmin > Vertex[Leaves[i].vertsLeaf[j]][2]) zmin = Vertex[Leaves[i].vertsLeaf[j]][2];
763                }
764        }
765
766        cx = (xmax + xmin)/2;
767        cy = (ymax + ymin)/2;
768        cz = (zmax + zmin)/2;
769
770        diameter = ((xmax-xmin)*(xmax-xmin)) + ((ymax-ymin)*(ymax-ymin)) + ((zmax-zmin)*(zmax-zmin));
771
772        return (diameter);
773}
774
775void TreeSimplifier::NormalizeDistance(float diameter)
776{
777        float dtmp;
778        for (int i=0; i<countLeaves;i++)
779        {
780                dtmp = Leaves[i].dist;
781                Leaves[i].dist = dtmp / diameter;
782               
783        }
784}
785/*
786
787        THIS FUNCTION IS OBSOLETE, AND NEEDS TO BE ERASED
788//--------------------------------------------------------------------------------------------------------------------------------
789//  establece el criterio como (K1 * dist + K2 * cop)/ K1+K2
790//--------------------------------------------------------------------------------------------------------------------------------
791void TreeSimplifier::EstableceCriterio(float diametro)
792{
793        //      Para empezar establezco K1 y K2 con 0,5
794        float coptmp2, coptmp, distmp2, distmp, criteriotmp;
795        int i, j, nhojasi, nhojasj;
796
797        int update_each = countLeaves / 30;
798       
799        for ( i=0; i<countLeaves;i++)
800        {
801                if (Leaves[i].existe == true)
802                {
803                        //incializo criterio a un numero elevado
804                        Leaves[i].criteria = 1000;
805                        nhojasi = int(Leaves[i].parentLeafCount);
806                        //coplanaridad
807                        for ( j =0; j<countLeaves;j++)
808                        {
809                                if (( j != i) && ( Leaves[j].existe == true)) // si la hoja aun existe
810                                {
811                                        //17/09/01 ANTES DE CALCULAR NADA, COMPRUEBO QUE ESTAS DOS HOJAS
812                                        // SE PODRIAN COLAPSAR, ED, QUE LA DIFERENCIA DE HOJAS QUE COLAPSAN
813                                        // ES COMO MÁXIMO 1
814
815                                        nhojasj = int(Leaves[j].parentLeafCount);
816
817                                        if ( abs((nhojasi - nhojasj)) < 2)
818                                        {
819                                                //coplanaridad y lo invierto
820                                                coptmp2 = Leaves[i].Coplanaridad(Leaves[j]);
821                                                coptmp = 1 - coptmp2;
822                                                //distancia y la normalizo
823                                                distmp2 = Hausdorff( Leaves[i], Leaves[j]);
824                                                distmp = distmp2 / diametro;
825                                                // calculo el criterio para esa hoja
826                                                criteriotmp = (( K1 * distmp * distmp ) + (K2 * coptmp * distmp))/ (K1 + K2);
827                                                //selecciono el criterio menor
828                                                if (Leaves[i].criteria > criteriotmp)
829                                                {
830                                                        Leaves[i].criteria = criteriotmp;
831                                                        Leaves[i].hoja_crit = j;
832                                                }
833                                        }
834                                }
835                               
836                        }
837                }
838        }
839}
840
841//--------------------------------------------------------------------------------------------------------------------------------
842//  establece el criterio despues de colapsar
843//--------------------------------------------------------------------------------------------------------------------------------
844void TreeSimplifier::EstableceCriterio2 ( float                 diametro,
845                                                                                                                                                                        long int        hojanueva)
846{
847        //Para empezar establezco K1 y K2 con 0,5
848        float  coptmp2, coptmp, distmp2, distmp, criteriotmp;
849        int i, j, nhojasi, nhojasj;
850        //ESTAN EN tres PROCEDIMIENTOS: COLAPSA, ESTABLECECRITERIO Y ESTABLECECRITERIO2
851       
852        for (i = 0; i < countLeaves; i++)
853        {
854                if ((Leaves[i].existe == true) && (i != hojanueva))
855                {       
856                        nhojasi = int(Leaves[i].parentLeafCount);
857                        //¿ SE HA DESACTIVADO LA HOJA_CRIT QUE GUARDABA LA HOJA?
858                        if ( Leaves[Leaves[i].hoja_crit].existe == false)
859                        {
860                                Leaves[i].criteria = 1000;
861
862                                //coplanaridad
863                                for ( j =0; j<countLeaves;j++)
864                                {
865                                        if (( j != i) && ( Leaves[j].existe == true)) // si la hoja aun existe
866                                        {
867                                                //17/09/01 ANTES DE CALCULAR NADA, COMPRUEBO QUE ESTAS DOS HOJAS
868                                                // SE PODRIAN COLAPSAR, ED, QUE LA DIFERENCIA DE HOJAS QUE COLAPSAN
869                                                // ES COMO MÁXIMO 1
870
871                                                nhojasj = int(Leaves[j].parentLeafCount);
872
873                                                if ( abs((nhojasi - nhojasj)) < 2)
874                                                {
875
876                                                        //coplanaridad y lo invierto
877                                                        coptmp2 = Leaves[i].Coplanaridad(Leaves[j]);
878                                                        coptmp = 1 - coptmp2;   
879                                                        //distancia y la normalizo     
880                                                //      distmp2 = Leaves[i].Distancia(Leaves[j]);
881                                                        distmp2 = Hausdorff( Leaves[i],  Leaves[j]);
882                                                        distmp = distmp2 / diametro;
883                                                        // calculo el criterio para esa hoja
884                                                        criteriotmp = (( K1 * distmp * distmp ) + (K2 * coptmp * distmp))/ (K1 + K2);
885                                                        //selecciono el criterio menor
886                                                        if (Leaves[i].criteria > criteriotmp)
887                                                        {
888                                                                Leaves[i].criteria = criteriotmp;
889                                                                Leaves[i].hoja_crit = j;
890                                                        }
891                                                }
892                                        }
893                                }
894                        }
895                        else
896                        { // CALCULARE SI EL CRITERIO CON ESTA HOJA ES MENOR QUE EL ANTERIOR
897                                nhojasj = int(Leaves[hojanueva].parentLeafCount);
898
899                                if ( abs((nhojasi - nhojasj)) < 2)
900                                {
901                                        //coplanaridad y lo invierto
902                                        coptmp2 = Leaves[i].Coplanaridad(Leaves[hojanueva]);
903                                        coptmp = 1 - coptmp2;
904                                        //distancia y la normalizo
905                                        //distmp2 = Leaves[i].Distancia(Leaves[hojanueva]);
906                                        distmp2 = Hausdorff( Leaves[i], Leaves[hojanueva]);
907                                        distmp = distmp2 / diametro;
908                                        // calculo el criterio para esa hoja
909                                        criteriotmp = (( K1 * distmp * distmp ) + (K2 * coptmp * distmp))/ (K1 + K2);
910                                        //selecciono el criterio menor
911                                        if (Leaves[i].criteria > criteriotmp)
912                                        {       
913                                                Leaves[i].criteria = criteriotmp;
914                                                Leaves[i].hoja_crit = hojanueva;
915                                        }
916                                }
917                        }
918                }
919        }
920}*/
921
922
923
924// Gets the laf with the minimum criteria
925long int TreeSimplifier::MinCriteria (void)
926{
927        float mincrit=0.0f;
928        long int which=-1;
929               
930/*      while (Leaves[i].existe != true) i++;
931        mincrit = Leaves[i].criteria;
932        which =i;
933        */
934        for (int i=0;i<countLeaves;i++)
935        {
936                if (Leaves[i].exists)
937                {
938                        if (mincrit>Leaves[i].criteria || which==-1)
939                        {
940                                mincrit = Leaves[i].criteria;
941                                which = i;
942                        }
943                }
944               
945        }
946        return which;
947}
948
949
950void TreeSimplifier::ChooseVertices(Leaf& leaf1, Leaf& leaf2, long int count)
951{
952        float a,b,c;
953        float dist[4];
954        int indices[4];
955       
956        a = Vertex[leaf1.vertsLeaf[0]][0];
957        b = Vertex[leaf1.vertsLeaf[0]][1];
958        c = Vertex[leaf1.vertsLeaf[0]][2];
959
960        dist[0] = ((leaf2.center[0]-a)*(leaf2.center[0]-a)) + ((leaf2.center[1]-b)*(leaf2.center[1]-b)) +
961                   ((leaf2.center[2]-c)*(leaf2.center[2]-c));
962
963       
964        a = Vertex[leaf1.vertsLeaf[1]][0];
965        b = Vertex[leaf1.vertsLeaf[1]][1];
966        c = Vertex[leaf1.vertsLeaf[1]][2];
967
968        dist[1] = ((leaf2.center[0]-a)*(leaf2.center[0]-a)) + ((leaf2.center[1]-b)*(leaf2.center[1]-b)) +
969                   ((leaf2.center[2]-c)*(leaf2.center[2]-c));
970
971        a = Vertex[leaf1.vertsLeaf[2]][0];
972        b = Vertex[leaf1.vertsLeaf[2]][1];
973        c = Vertex[leaf1.vertsLeaf[2]][2];
974
975        dist[2] = ((leaf2.center[0]-a)*(leaf2.center[0]-a)) + ((leaf2.center[1]-b)*(leaf2.center[1]-b)) +
976                   ((leaf2.center[2]-c)*(leaf2.center[2]-c));
977
978
979        a = Vertex[leaf1.vertsLeaf[3]][0];
980        b = Vertex[leaf1.vertsLeaf[3]][1];
981        c = Vertex[leaf1.vertsLeaf[3]][2];
982
983        dist[3] = ((leaf2.center[0]-a)*(leaf2.center[0]-a)) + ((leaf2.center[1]-b)*(leaf2.center[1]-b)) +
984                   ((leaf2.center[2]-c)*(leaf2.center[2]-c));
985
986        for ( int i=0;i<4;i++) indices[i]=i;
987
988        TwoGreater(dist, indices);
989       
990        Leaves[countLeaves].vertsLeaf[0] = leaf1.vertsLeaf[indices[0]];
991        Leaves[countLeaves].vertsLeaf[1] = leaf1.vertsLeaf[indices[1]];
992       
993
994       
995        a = Vertex[leaf2.vertsLeaf[0]][0];
996        b = Vertex[leaf2.vertsLeaf[0]][1];
997        c = Vertex[leaf2.vertsLeaf[0]][2];
998
999        dist[0] = ((leaf1.center[0]-a)*(leaf1.center[0]-a)) + ((leaf1.center[1]-b)*(leaf1.center[1]-b)) +
1000                   ((leaf1.center[2]-c)*(leaf1.center[2]-c));
1001
1002       
1003        a = Vertex[leaf2.vertsLeaf[1]][0];
1004        b = Vertex[leaf2.vertsLeaf[1]][1];
1005        c = Vertex[leaf2.vertsLeaf[1]][2];
1006
1007        dist[1] = ((leaf2.center[0]-a)*(leaf2.center[0]-a)) + ((leaf1.center[1]-b)*(leaf1.center[1]-b)) +
1008                   ((leaf2.center[2]-c)*(leaf2.center[2]-c));
1009
1010        a = Vertex[leaf2.vertsLeaf[2]][0];
1011        b = Vertex[leaf2.vertsLeaf[2]][1];
1012        c = Vertex[leaf2.vertsLeaf[2]][2];
1013
1014        dist[2] = ((leaf1.center[0]-a)*(leaf1.center[0]-a)) + ((leaf1.center[1]-b)*(leaf1.center[1]-b)) +
1015                   ((leaf1.center[2]-c)*(leaf1.center[2]-c));
1016
1017
1018        a = Vertex[leaf2.vertsLeaf[3]][0];
1019        b = Vertex[leaf2.vertsLeaf[3]][1];
1020        c = Vertex[leaf2.vertsLeaf[3]][2];
1021
1022        dist[3] = ((leaf1.center[0]-a)*(leaf1.center[0]-a)) + ((leaf1.center[1]-b)*(leaf1.center[1]-b)) +
1023                   ((leaf1.center[2]-c)*(leaf1.center[2]-c));
1024
1025        for ( int i=0;i<4;i++) indices[i]=i;
1026
1027        TwoGreater(dist, indices);
1028       
1029        Leaves[countLeaves].vertsLeaf[2] = leaf2.vertsLeaf[indices[0]];
1030        Leaves[countLeaves].vertsLeaf[3] = leaf2.vertsLeaf[indices[1]];
1031
1032
1033
1034       
1035}
1036
1037// Add leaves
1038long int TreeSimplifier::Collapse(float diam)
1039{
1040        long int i=0, which=-1;
1041        long int other = -1;
1042        float coptmp, coptmp2, distmp, distmp2, criteriatmp;
1043
1044        which=MinCriteria();
1045        other = Leaves[which].leafCrit;
1046
1047        //desactivo las hojas cercanas
1048        Leaves[which].exists = false;
1049        Leaves[other].exists = false;
1050        //creo la hoja nueva
1051
1052        Leaves[countLeaves].leafNear = -1;
1053        Leaves[countLeaves].dist = 0;
1054        Leaves[countLeaves].leafCop = -1;
1055        Leaves[countLeaves].coplanar = 0;
1056
1057        Leaves[countLeaves].parentLeafCount = Leaves[which].parentLeafCount + Leaves[other].parentLeafCount;
1058//      Leaves[countLeaves].parentLeafCount = 10;
1059        ChooseVertices(Leaves[which], Leaves[other], countLeaves);
1060        CalculateLeafCenter (Leaves[countLeaves]);
1061        CalculateLeafNormal (Leaves[countLeaves]);
1062
1063        // find out the minimal octree node that can contain this leaf
1064        LeafOctree *onode = GetMinOctreeNodeForLeaf(octree,Leaves[countLeaves]);
1065        octree_owning_leaf[countLeaves] = onode;
1066        onode->leaves.push_back(countLeaves);
1067
1068/*      if (Leaves[countLeaves].parentLeafCount > 60 )
1069        {
1070                        Leaves[countLeaves].existe = false;
1071                        activeLeaves--;
1072        }
1073        else
1074        {               */
1075        Leaves[countLeaves].exists = true;
1076                Leaves[countLeaves].criteria = 1000000.0f;
1077                Leaves[countLeaves].idTriangle[0] = countLeaves*2;
1078                Leaves[countLeaves].idTriangle[1] = countLeaves*2+1;
1079
1080                float area_leaf_i = CalculateLeafArea(Leaves[i])/diam;
1081
1082                for (int j = 0; j < (countLeaves+1); j++){
1083/*              int visit_parents = VISIT_PARENTS_DEEP;
1084                for (LeafOctree *octreenode = octree_owning_leaf[i]; octreenode && visit_parents>=0; octreenode=octreenode->parent, visit_parents--)
1085                {
1086                        for (std::vector<int>::iterator it=octreenode->leaves.begin(); it!=octreenode->leaves.end(); it++)
1087                        {
1088                                int j = *it;*/
1089
1090                                if (j != countLeaves && Leaves[j].exists)
1091                                {
1092                                        coptmp2 = Leaves[countLeaves].Coplanarity(Leaves[j]);
1093                                        coptmp = 1 - coptmp2;
1094                                        //distmp2 = HausdorffOptimized(Leaves[countLeaves], Leaves[j]);
1095                                        distmp2 = DistanceFromCenters(Leaves[countLeaves], Leaves[j]);
1096                                        distmp = distmp2 / diam;
1097                                       
1098                                        criteriatmp = (( K1 * distmp * distmp ) + (K2 * coptmp * distmp))/ (K1 + K2);
1099                                        criteriatmp *= CalculateLeafArea(Leaves[j])/diam + area_leaf_i;
1100                                        criteriatmp *= Leaves[j].parentLeafCount + Leaves[i].parentLeafCount;
1101                                       
1102                                        if (Leaves[countLeaves].criteria > criteriatmp)
1103                                        {
1104                                                Leaves[countLeaves].criteria = criteriatmp;
1105                                                Leaves[countLeaves].leafCrit = j;
1106                                        }
1107                                }
1108        //              }
1109                }
1110//      }       
1111
1112        // Crear el paso de simplificación
1113        Geometry::TreeSimplificationSequence::Step pasosimp;
1114        pasosimp.mV0=Leaves[which].idTriangle[0];
1115        pasosimp.mV1=Leaves[which].idTriangle[1];
1116        pasosimp.mT0=Leaves[other].idTriangle[0];
1117        pasosimp.mT1=Leaves[other].idTriangle[1];
1118
1119        // Nuevos vértices
1120        pasosimp.mNewQuad[0]=Leaves[countLeaves].vertsLeaf[0];
1121        pasosimp.mNewQuad[1]=Leaves[countLeaves].vertsLeaf[1];
1122        pasosimp.mNewQuad[2]=Leaves[countLeaves].vertsLeaf[2];
1123        pasosimp.mNewQuad[3]=Leaves[countLeaves].vertsLeaf[3];
1124
1125        // Insertar el paso de simplificación
1126        mtreesimpsequence->mSteps.push_back(pasosimp);
1127
1128        //incremento el numero de hojas
1129        countLeaves++;
1130        activeLeaves--;
1131        return (countLeaves-1);//porque lo he incrementado en la linea de antes
1132}
1133
1134void TreeSimplifier::BuildOutputMesh(int idMeshLeaves)
1135{
1136        // Calculate the resulting index count after simplifying
1137        long int tamIndex       =       activeLeaves * 6; // Each leaf has got 6 indices
1138
1139        delete [] mesh->mSubMesh[idMeshLeaves].mIndex;
1140        mesh->mSubMesh[idMeshLeaves].mIndexCount=tamIndex;
1141        mesh->mSubMesh[idMeshLeaves].mIndex=new Geometry::Index[tamIndex];
1142
1143        // Iterate through all active leaves and dump them to the final mesh
1144        Geometry::Index index=0;
1145
1146        int update_each = countLeaves / 10;
1147       
1148        for (long j=0; j<countLeaves;j++)
1149        {
1150                static int ticks_since_last_update = 0;
1151                if (mUPB && ticks_since_last_update>=update_each)
1152                {
1153                        ticks_since_last_update=0;
1154                        mUPB(1);
1155                }
1156                ticks_since_last_update++;
1157               
1158                if (Leaves[j].exists)
1159                {
1160//                      if (index<6)
1161//                      {
1162                                mesh->mSubMesh[idMeshLeaves].mIndex[index]=Leaves[j].vertsLeaf[0];
1163                                mesh->mSubMesh[idMeshLeaves].mIndex[index+1]=Leaves[j].vertsLeaf[1];
1164                                mesh->mSubMesh[idMeshLeaves].mIndex[index+2]=Leaves[j].vertsLeaf[2];
1165                                mesh->mSubMesh[idMeshLeaves].mIndex[index+3]=Leaves[j].vertsLeaf[2];
1166                                mesh->mSubMesh[idMeshLeaves].mIndex[index+4]=Leaves[j].vertsLeaf[1];
1167                                mesh->mSubMesh[idMeshLeaves].mIndex[index+5]=Leaves[j].vertsLeaf[3];
1168/*                      }
1169                        else
1170                        {
1171                                mesh->mSubMesh[idMeshLeaves].mIndex[index]=0;
1172                                mesh->mSubMesh[idMeshLeaves].mIndex[index+1]=0;
1173                                mesh->mSubMesh[idMeshLeaves].mIndex[index+2]=0;
1174                                mesh->mSubMesh[idMeshLeaves].mIndex[index+3]=0;
1175                                mesh->mSubMesh[idMeshLeaves].mIndex[index+4]=0;
1176                                mesh->mSubMesh[idMeshLeaves].mIndex[index+5]=0;
1177                        }*/
1178                        index=index+6;
1179                }
1180        }
1181}
1182
1183LeafOctree* TreeSimplifier::CreateLeafOctree(int deep)
1184{
1185        LeafOctree *leafoctree = new LeafOctree;
1186        leafoctree->parent=NULL;
1187        leafoctree->left = Vertex[0][0];
1188        leafoctree->right = Vertex[0][0];
1189        leafoctree->bottom = Vertex[0][1];
1190        leafoctree->top = Vertex[0][1];
1191        leafoctree->front = Vertex[0][2];
1192        leafoctree->back = Vertex[0][2];
1193
1194        // calcualte the limits of the octree
1195        for (size_t     i       =       0; i < vertex_count; i++)
1196        {
1197                if (leafoctree->left>Vertex[i][0])
1198                        leafoctree->left=Vertex[i][0];
1199                if (leafoctree->right<Vertex[i][0])
1200                        leafoctree->right=Vertex[i][0];
1201                if (leafoctree->bottom>Vertex[i][1])
1202                        leafoctree->bottom=Vertex[i][1];
1203                if (leafoctree->top<Vertex[i][1])
1204                        leafoctree->top=Vertex[i][1];
1205                if (leafoctree->front>Vertex[i][2])
1206                        leafoctree->front=Vertex[i][2];
1207                if (leafoctree->back<Vertex[i][2])
1208                        leafoctree->back=Vertex[i][2];
1209        }
1210
1211        // setup the initial leaves buffer
1212        leafoctree->leaves.clear();
1213        octree_owning_leaf=new LeafOctree*[countLeaves*2]; // *2 to reserve memory for all simplified leaves
1214        for (int i=0; i<countLeaves; i++)
1215        {
1216                leafoctree->leaves.push_back(i);
1217                octree_owning_leaf[i] = leafoctree;
1218        }
1219
1220        // generate the octree given an arbitrary depth
1221        RecursiveCreateLeafOctree(leafoctree,deep);
1222
1223        // fill the octree
1224        RecursiveFillOctreeWithLeaves(leafoctree);
1225
1226        return leafoctree;
1227}
1228
1229void TreeSimplifier::RecursiveCreateLeafOctree(LeafOctree *parent, int deep)
1230{
1231        for (int i=0; i<8; i++)
1232        {
1233                parent->children[i] = NULL;
1234                if (deep>0)
1235                {
1236                        LeafOctree *child = parent->children[i] = new LeafOctree;
1237
1238                        child->parent = parent;
1239                        switch(i)
1240                        {
1241                        case 0:
1242                                child->left=parent->left;
1243                                child->right=(parent->left+parent->right)*0.5f;
1244                                child->bottom=(parent->bottom+parent->top)*0.5f;
1245                                child->top=parent->top;
1246                                child->front=parent->front;
1247                                child->back=(parent->front+parent->back)*0.5f; break;
1248                        case 1:
1249                                child->left=(parent->left+parent->right)*0.5f;
1250                                child->right=parent->right;
1251                                child->bottom=(parent->bottom+parent->top)*0.5f;
1252                                child->top=parent->top;
1253                                child->front=parent->front;
1254                                child->back=(parent->front+parent->back)*0.5f; break;
1255                        case 2:
1256                                child->left=parent->left;
1257                                child->right=(parent->left+parent->right)*0.5f;
1258                                child->bottom=parent->bottom;
1259                                child->top=(parent->bottom+parent->top)*0.5f;
1260                                child->front=parent->front;
1261                                child->back=(parent->front+parent->back)*0.5f; break;
1262                        case 3:
1263                                child->left=(parent->left+parent->right)*0.5f;
1264                                child->right=parent->right;
1265                                child->bottom=parent->bottom;
1266                                child->top=(parent->bottom+parent->top)*0.5f;
1267                                child->front=parent->front;
1268                                child->back=(parent->front+parent->back)*0.5f; break;
1269                        case 4:
1270                                child->left=parent->left;
1271                                child->right=(parent->left+parent->right)*0.5f;
1272                                child->bottom=(parent->bottom+parent->top)*0.5f;
1273                                child->top=parent->top;
1274                                child->back=parent->back;
1275                                child->front=(parent->front+parent->back)*0.5f; break;
1276                        case 5:
1277                                child->left=(parent->left+parent->right)*0.5f;
1278                                child->right=parent->right;
1279                                child->bottom=(parent->bottom+parent->top)*0.5f;
1280                                child->top=parent->top;
1281                                child->back=parent->back;
1282                                child->front=(parent->front+parent->back)*0.5f; break;
1283                        case 6:
1284                                child->left=parent->left;
1285                                child->right=(parent->left+parent->right)*0.5f;
1286                                child->bottom=parent->bottom;
1287                                child->top=(parent->bottom+parent->top)*0.5f;
1288                                child->back=parent->back;
1289                                child->front=(parent->front+parent->back)*0.5f; break;
1290                        case 7:
1291                                child->left=(parent->left+parent->right)*0.5f;
1292                                child->right=parent->right;
1293                                child->bottom=parent->bottom;
1294                                child->top=(parent->bottom+parent->top)*0.5f;
1295                                child->back=parent->back;
1296                                child->front=(parent->front+parent->back)*0.5f; break;
1297                        }
1298                        RecursiveCreateLeafOctree(parent->children[i],deep-1);
1299                }
1300        }
1301
1302}
1303
1304void TreeSimplifier::RecursiveFillOctreeWithLeaves(LeafOctree *o)
1305{   
1306        for (int i=0; i<8; i++)
1307        {
1308                LeafOctree *child = o->children[i];
1309                if (child)
1310                {
1311                        for (std::vector<int>::iterator it=o->leaves.begin(); it!=o->leaves.end(); )
1312                        {
1313                                int idleaf = *it;
1314                                int idv0 = Leaves[idleaf].vertsLeaf[0];
1315                                int idv1 = Leaves[idleaf].vertsLeaf[1];
1316                                int idv2 = Leaves[idleaf].vertsLeaf[2];
1317                                int idv3 = Leaves[idleaf].vertsLeaf[3];
1318
1319                                float * v0 = Vertex[idv0];
1320                                float * v1 = Vertex[idv1];
1321                                float * v2 = Vertex[idv2];
1322                                float * v3 = Vertex[idv3];
1323
1324                                // if the leaf fits completely inside a child, the child owns it
1325                                if (child->PointInsideOctree(v0[0],v0[1],v0[2]) &&
1326                                        child->PointInsideOctree(v1[0],v1[1],v1[2]) &&
1327                                        child->PointInsideOctree(v2[0],v2[1],v2[2]) &&
1328                                        child->PointInsideOctree(v3[0],v3[1],v3[2]))
1329                                {
1330                                        child->leaves.push_back(idleaf);
1331                                        it = o->leaves.erase(it);
1332                                        octree_owning_leaf[idleaf] = child;
1333                                }
1334                                else
1335                                        it++;
1336                        }
1337                        // recurse all valid children
1338                        RecursiveFillOctreeWithLeaves(o->children[i]);
1339                }
1340        }
1341}
1342
1343void TreeSimplifier::SetCriteriaOptimized(float diam)
1344{
1345        float coptmp2, coptmp, distmp2, distmp, criteriatmp;
1346        int i, j, nleavesi, nleavesj;
1347
1348        int update_each = countLeaves / 20;
1349       
1350        for ( i=0; i<countLeaves;i++)
1351        {
1352                static int ticks_since_last_update = 0;
1353                if (mUPB && ticks_since_last_update>=update_each)
1354                {
1355                        ticks_since_last_update=0;
1356                        mUPB(1);
1357                }
1358                ticks_since_last_update++;
1359               
1360                if (Leaves[i].exists == true)
1361                {
1362                        Leaves[i].criteria = 1000000.0f;
1363                        nleavesi = int(Leaves[i].parentLeafCount);
1364                        float area_leaf_i = CalculateLeafArea(Leaves[i])/diam;
1365
1366                        // Use only the leaves which belong the the same octree node as leaf i or to its parents               
1367                        int visit_parents = VISIT_PARENTS_DEEP;
1368                        for (LeafOctree *octreenode = octree_owning_leaf[i]; octreenode; octreenode=octreenode->parent, visit_parents--)
1369                        {
1370                                for (std::vector<int>::iterator it=octreenode->leaves.begin(); it!=octreenode->leaves.end(); it++)
1371                                {
1372                                        j = *it;
1373//                              for (j=0; j<countLeaves; j++)
1374//                              {
1375                                        if (j!=i && Leaves[j].exists)
1376                                        {
1377                                                nleavesj = int(Leaves[j].parentLeafCount);
1378
1379        //                                      if ( abs((nleavesi - nleavesj)) < 2)
1380        //                                      {                                                       
1381                                                        coptmp2 = Leaves[i].Coplanarity(Leaves[j]);
1382                                                        coptmp = 1 - coptmp2;
1383
1384                                                        //distmp2 = HausdorffOptimized( Leaves[i], Leaves[j]);
1385                                                        distmp2 = DistanceFromCenters(Leaves[i], Leaves[j]);
1386                                                        distmp = distmp2 / diam;
1387
1388                                                        //criteriatmp = (( K1 * distmp * distmp ) + (K2 * coptmp * distmp))/ (K1 + K2);
1389                                                        criteriatmp = distmp;
1390                                                        // select the lowest
1391                                                        criteriatmp *= CalculateLeafArea(Leaves[j])/diam + area_leaf_i;
1392                                                        criteriatmp *= Leaves[j].parentLeafCount + Leaves[i].parentLeafCount;
1393                                                        if (Leaves[i].criteria > criteriatmp)
1394                                                        {
1395                                                                Leaves[i].criteria = criteriatmp;
1396                                                                Leaves[i].leafCrit = j;
1397                                                        }
1398                //                              }
1399                                        }
1400                                }                               
1401                        }
1402                }
1403        }
1404}
1405
1406
1407void TreeSimplifier::SetCriteria2Optimized(float diam, long int newleaf)
1408{
1409        float  coptmp2, coptmp, distmp2, distmp, criteriatmp;
1410        int i, j, nleavesi, nleavesj;
1411       
1412        for (i = 0; i < countLeaves; i++)
1413        {
1414                if (Leaves[i].exists && i!=newleaf)
1415                {       
1416                        float area_leaf_i = CalculateLeafArea(Leaves[i])/diam;
1417                        nleavesi = int(Leaves[i].parentLeafCount);
1418                        if ( Leaves[Leaves[i].leafCrit].exists == false)
1419                        {
1420                                Leaves[i].criteria = 1000000.0f;
1421
1422                                int visit_parents = VISIT_PARENTS_DEEP;
1423                                for (LeafOctree *octreenode = octree_owning_leaf[i]; octreenode; octreenode=octreenode->parent, visit_parents--)
1424                                {
1425                                        for (std::vector<int>::iterator it=octreenode->leaves.begin(); it!=octreenode->leaves.end(); it++)
1426                                        {
1427                                                j = *it;
1428        //                      for (j=0; j<countLeaves; j++)
1429        //                      {
1430                                                if (j!=i && Leaves[j].exists)
1431                                                {
1432                                                        nleavesj = int(Leaves[j].parentLeafCount);
1433
1434//                                                      if ( abs((nleavesi - nleavesj)) < 2)
1435//                                                      {
1436                                                                coptmp2 = Leaves[i].Coplanarity(Leaves[j]);
1437                                                                coptmp = 1 - coptmp2;   
1438                                                       
1439                                                       
1440                                                                //distmp2 = HausdorffOptimized( Leaves[i],  Leaves[j]);
1441                                                                distmp2 = DistanceFromCenters( Leaves[i], Leaves[j]);
1442                                                                distmp = distmp2 / diam;
1443
1444                                                                criteriatmp = (( K1 * distmp * distmp ) + (K2 * coptmp * distmp))/ (K1 + K2);
1445                                                                criteriatmp *= CalculateLeafArea(Leaves[j])/diam + area_leaf_i;
1446                                                                criteriatmp *= Leaves[j].parentLeafCount + Leaves[i].parentLeafCount;
1447                                                               
1448                                                                // select the leaf with the lowest criteria
1449                                                                if (Leaves[i].criteria > criteriatmp)
1450                                                                {
1451                                                                        Leaves[i].criteria = criteriatmp;
1452                                                                        Leaves[i].leafCrit = j;
1453                                                                }
1454                                                //      }
1455                                                }
1456                                        }
1457                                }
1458                        }
1459                        else
1460                        {
1461                                nleavesj = int(Leaves[newleaf].parentLeafCount);
1462
1463//                              if ( abs((nleavesi - nleavesj)) < 2)
1464//                              {
1465                                        coptmp2 = Leaves[i].Coplanarity(Leaves[newleaf]);
1466                                        coptmp = 1 - coptmp2;
1467
1468                                        //distmp2 = HausdorffOptimized( Leaves[i], Leaves[newleaf]);
1469                                        distmp2 = DistanceFromCenters( Leaves[i], Leaves[newleaf]);
1470                                        distmp = distmp2 / diam;
1471                                       
1472                                        criteriatmp = (( K1 * distmp * distmp ) + (K2 * coptmp * distmp))/ (K1 + K2);
1473                                        criteriatmp *= CalculateLeafArea(Leaves[newleaf])/diam + area_leaf_i;
1474                                        criteriatmp *= Leaves[newleaf].parentLeafCount + Leaves[i].parentLeafCount;
1475                                       
1476                                        if (Leaves[i].criteria > criteriatmp)
1477                                        {       
1478                                                Leaves[i].criteria = criteriatmp;
1479                                                Leaves[i].leafCrit = newleaf;
1480                                        }
1481//                              }
1482                        }
1483                }
1484        }
1485}
1486
1487
1488LeafOctree *TreeSimplifier::GetMinOctreeNodeForLeaf(LeafOctree *start, const Leaf &leaf)
1489{
1490        int idv0 = leaf.vertsLeaf[0];
1491        int idv1 = leaf.vertsLeaf[1];
1492        int idv2 = leaf.vertsLeaf[2];
1493        int idv3 = leaf.vertsLeaf[3];
1494
1495        float * v0 = Vertex[idv0];
1496        float * v1 = Vertex[idv1];
1497        float * v2 = Vertex[idv2];
1498        float * v3 = Vertex[idv3];
1499
1500        if (start->PointInsideOctree(v0[0],v0[1],v0[2]) &&
1501                start->PointInsideOctree(v1[0],v1[1],v1[2]) &&
1502                start->PointInsideOctree(v2[0],v2[1],v2[2]) &&
1503                start->PointInsideOctree(v3[0],v3[1],v3[2]))
1504        {
1505                return start;
1506        }
1507
1508        for (int i=0; i<8; i++)
1509                if (start->children[i])
1510                {
1511                        LeafOctree *selectedNode = GetMinOctreeNodeForLeaf(start->children[i],leaf);
1512                        if (selectedNode)
1513                                return selectedNode;
1514                }
1515
1516        return NULL;
1517}
1518
1519float TreeSimplifier::DistanceFromCenters(const Leaf &leaf1, const Leaf &leaf2) const
1520{
1521        float onex, oney, onez;
1522        float twox, twoy, twoz;
1523        float threex, threey, threez;
1524        float fourx, foury, fourz;
1525        float x1, y1, z1;
1526        float x2, y2, z2;
1527        float x3, y3, z3;
1528        float x4, y4, z4;
1529
1530        onex   = Vertex[leaf1.vertsLeaf[0]][0]; oney   = Vertex[leaf1.vertsLeaf[0]][1]; onez   = Vertex[leaf1.vertsLeaf[0]][2];
1531        twox   = Vertex[leaf1.vertsLeaf[1]][0]; twoy   = Vertex[leaf1.vertsLeaf[1]][1]; twoz   = Vertex[leaf1.vertsLeaf[1]][2];
1532        threex = Vertex[leaf1.vertsLeaf[2]][0]; threey = Vertex[leaf1.vertsLeaf[2]][1]; threez = Vertex[leaf1.vertsLeaf[2]][2];
1533        fourx  = Vertex[leaf1.vertsLeaf[3]][0]; foury  = Vertex[leaf1.vertsLeaf[3]][1]; fourz  = Vertex[leaf1.vertsLeaf[3]][2];
1534
1535        float center1x = (onex + twox + threex + fourx)*0.25f;
1536        float center1y = (oney + twoy + threey + foury)*0.25f;
1537        float center1z = (onez + twoz + threez + fourz)*0.25f;
1538
1539        x1 = Vertex[leaf2.vertsLeaf[0]][0]; y1 = Vertex[leaf2.vertsLeaf[0]][1]; z1 = Vertex[leaf2.vertsLeaf[0]][2];
1540        x2 = Vertex[leaf2.vertsLeaf[1]][0]; y2 = Vertex[leaf2.vertsLeaf[1]][1]; z2 = Vertex[leaf2.vertsLeaf[1]][2];
1541        x3 = Vertex[leaf2.vertsLeaf[2]][0]; y3 = Vertex[leaf2.vertsLeaf[2]][1]; z3 = Vertex[leaf2.vertsLeaf[2]][2];
1542        x4 = Vertex[leaf2.vertsLeaf[3]][0]; y4 = Vertex[leaf2.vertsLeaf[3]][1]; z4 = Vertex[leaf2.vertsLeaf[3]][2];
1543
1544        float center2x = (x1 + x2 + x3 + x4)*0.25f;
1545        float center2y = (y1 + y2 + y3 + y4)*0.25f;
1546        float center2z = (z1 + z2 + z3 + z4)*0.25f;
1547
1548        return distan(center1x,center1y,center1z,center2x,center2y,center2z);
1549}
1550
1551
1552bool TreeSimplifier::PruneOctree(LeafOctree *octree)
1553{
1554        if (octree->HasChildren())
1555        {
1556                for (int i=0; i<8; i++)
1557                        if (octree->children[i])
1558                                if (PruneOctree(octree->children[i])) // if child was pruned, nullify it
1559                                        octree->children[i]=NULL;                       
1560        }
1561        else
1562        {
1563                if (octree->parent)
1564                {
1565                        for (std::vector<int>::iterator it=octree->leaves.begin(); it!=octree->leaves.end(); it++)
1566                        {
1567                                int idleaf = *it;
1568                                octree->parent->leaves.push_back(idleaf);
1569                                octree_owning_leaf[idleaf] = octree->parent;
1570                        }
1571                        return true;
1572                }
1573        }
1574        return false;
1575}
1576
1577void TreeSimplifier::CrossProduct(const float *v1, const float *v2, float *res) const
1578{
1579        res[0] = v1[1]*v2[2] - v1[2]*v2[1];
1580        res[1] = v1[2]*v2[0] - v1[0]*v2[2];
1581        res[2] = v1[0]*v2[1] - v1[1]*v2[0];
1582}
1583
1584float TreeSimplifier::SquaredModule(const float *v) const
1585{
1586        return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
1587}
1588
1589
1590float TreeSimplifier::CalculateLeafArea(Leaf &leaf) const
1591{
1592        int idv0 = leaf.vertsLeaf[0];
1593        int idv1 = leaf.vertsLeaf[1];
1594        int idv2 = leaf.vertsLeaf[2];
1595        int idv3 = leaf.vertsLeaf[3];
1596
1597        float * v0 = Vertex[idv0];
1598        float * v1 = Vertex[idv1];
1599        float * v2 = Vertex[idv2];
1600        float * v3 = Vertex[idv3];
1601
1602        float v1v0[3] = {v1[0]-v0[0],v1[1]-v0[1],v1[2]-v0[2]};
1603        float v2v0[3] = {v2[0]-v0[0],v2[1]-v0[1],v2[2]-v0[2]};
1604        float v2v1[3] = {v2[0]-v1[0],v2[1]-v1[1],v2[2]-v1[2]};
1605        float v3v1[3] = {v3[0]-v1[0],v3[1]-v1[1],v3[2]-v1[2]};
1606       
1607        float cross1[3], cross2[3];
1608
1609        CrossProduct(v1v0,v2v0,cross1);
1610        CrossProduct(v2v1,v3v1,cross2);
1611
1612        float mod1 = SquaredModule(cross1);
1613        float mod2 = SquaredModule(cross2);
1614
1615        return mod1*0.5f + mod2*0.5f;
1616}
Note: See TracBrowser for help on using the repository browser.