source: GTP/trunk/App/Demos/Illum/PathMap/Wsimplex.cpp @ 896

Revision 896, 12.3 KB checked in by szirmay, 19 years ago (diff)
RevLine 
[896]1// Wsimplex.cpp: implementation of the Wsimplex class.
2//
3//////////////////////////////////////////////////////////////////////
4
5#include "dxstdafx.h"
6#include "Wsimplex.h"
7#include "Vector.hpp"
8#include "Wbucket.h"
9#include "WsubMesh.h"
10
11//////////////////////////////////////////////////////////////////////
12// ruction/Destruction
13//////////////////////////////////////////////////////////////////////
14
15Wsimplex::Wsimplex(int dim)
16{
17        deleted = 0;
18        circumcentre = new Vector(0,0,0);
19        dimension = dim;
20
21        vertex = new Vector*[dimension+1];
22        neighbour = new Wsimplex*[dimension+1];
23}
24
25Wsimplex::Wsimplex(Vector *q, Vector *w, int order)
26{       
27        dimension = 1;
28        deleted = 0;
29        vertex = new Vector*[dimension+1];
30        neighbour = new Wsimplex*[dimension+1];
31
32        circumcentre = new Vector(0,0,0);
33
34        neighbour[0] = NULL;
35        neighbour[1] = NULL;
36        inFrontOf = NULL;
37        behind = NULL;
38
39        if(!order)
40        {
41        vertex[0] = q;
42        vertex[1] = w;
43        }
44        else
45        {
46                if(q<w)
47                {
48                        vertex[0]=q;
49                        vertex[1]=w;
50                }
51                else
52                {
53                        vertex[0] = w;
54                        vertex[1] = q;
55                }
56
57        }
58//WARNING correct ordering by pointer value assumed if(order == 0)
59        updateCentre();
60}
61
62
63Wsimplex::Wsimplex(Vector *q, Vector *w, Vector *e, int order)
64{       
65        dimension = 2;
66        deleted = 0;
67        vertex = new Vector*[dimension+1];
68        neighbour = new Wsimplex*[dimension+1];
69
70        circumcentre = new Vector(0,0,0);
71        neighbour[0] = NULL;
72        neighbour[1] = NULL;
73        neighbour[2] = NULL;
74        inFrontOf = NULL;
75        behind = NULL;
76
77        if(!order)
78        {
79        vertex[0] = q;
80        vertex[1] = w;
81        vertex[2] = e;
82        }
83        else
84        {
85                if(q<w)
86                {
87                        if(q<e)
88                        {
89                                vertex[0]=q;
90                                if(w<e)
91                                {
92                                        vertex[1]=w;
93                                        vertex[2]=e;
94                                }
95                                else
96                                {
97                                        vertex[1]=e;
98                                        vertex[2]=w;
99                                }
100                        }
101                        else
102                        {
103                                vertex[0] = e;
104                                vertex[1] = q;
105                                vertex[2] = w;
106                        }
107                }
108                else
109                {
110                        if(q<e)
111                        {
112                                vertex[0] = w;
113                                vertex[1] = q;
114                                vertex[2] = e;
115                        }
116                        else
117                        {
118                                vertex[2] = q;
119                                if(w<e)
120                                {
121                                        vertex[0]=w;
122                                        vertex[1]=e;
123                                }
124                                else
125                                {
126                                        vertex[0]=e;
127                                        vertex[1]=w;
128                                }
129                        }
130                }
131        }
132
133//WARNING correct ordering by pointer value assumed if(order == 0)
134        updateCentre();
135}
136
137Wsimplex::Wsimplex(Vector *sharp, Wsimplex *flat)
138{
139        dimension = flat->getDimension()+1;
140        if(dimension == 2)
141                printf("brek \n");
142
143        vertex = new Vector*[dimension+1];
144        neighbour = new Wsimplex*[dimension+1];
145       
146        deleted = 0;
147        circumcentre = new Vector(0,0,0);
148       
149        for(int shg = 0;shg < dimension+1;shg++)
150                neighbour[shg] = NULL;
151        inFrontOf = NULL;
152        behind = NULL;
153
154        Vector** tri = flat->getVertices();
155
156        int sindex;
157        for(sindex = 0;sindex<dimension;sindex++)
158                if(sharp < tri[sindex]) break;;
159
160        vertex[sindex] = sharp;
161        neighbour[sindex] = flat->getBehind();
162
163        switch(dimension)
164        {
165        case 3: vertex[(sindex>2)?2:3] = tri[2];
166                        neighbour[(sindex>2)?2:3] = flat->getNeighboursFront(2);
167        case 2: vertex[(sindex>1)?1:2] = tri[1];
168                        neighbour[(sindex>1)?1:2] = flat->getNeighboursFront(1);
169        case 1: vertex[sindex?0:1] = tri[0];
170                        neighbour[sindex?0:1] = flat->getNeighboursFront(0);
171        }
172       
173//tell the neighbours to link here
174        for(int i=0;i<dimension+1;i++)
175                if(neighbour[i])
176                        neighbour[i]->helloNeighbour(this);
177
178        updateCentre();
179}
180
181
182Wsimplex::~Wsimplex()
183{
184        if(!backlist.empty())
185        {
186                ::std::vector<Wbucket*>::const_iterator a = backlist.begin();
187                while(a != backlist.end())
188                {
189                        (*a)->eliminate(this);
190                        a++;
191                }
192        }
193//WARNING cc is possibly used in refine
194        delete circumcentre;
195}
196
197
198Vector* Wsimplex::getCircumcentre()
199{
200        return circumcentre;
201}
202
203double Wsimplex::getRadius()
204{
205        return radius;
206}
207
208int Wsimplex::inCircle(Vector *guess)
209{
210        Vector dist = *circumcentre - *guess;
211        return sqrad - 0.0000001> dist.norm2();
212}
213
214void Wsimplex::linkback(Wbucket *container)
215{
216        backlist.push_back(container);
217}
218
219
220void Wsimplex::updateCentre()
221{
222        if(dimension == 3)
223        {
224        double bax = vertex[1]->x - vertex[0]->x;
225        double cax = vertex[2]->x - vertex[0]->x;
226        double dax = vertex[3]->x - vertex[0]->x;
227        double bay = vertex[1]->y - vertex[0]->y;
228        double cay = vertex[2]->y - vertex[0]->y;
229        double day = vertex[3]->y - vertex[0]->y;
230        double baz = vertex[1]->z - vertex[0]->z;
231        double caz = vertex[2]->z - vertex[0]->z;
232        double daz = vertex[3]->z - vertex[0]->z;
233
234        double bxcx = bay * caz - baz * cay;
235        double bxcy = baz * cax - bax * caz;
236        double bxcz = bax * cay - bay * cax;
237
238        double dxbx = day * baz - daz * bay;
239        double dxby = daz * bax - dax * baz;
240        double dxbz = dax * bay - day * bax;
241
242        double cxdx = cay * daz - caz * day;
243        double cxdy = caz * dax - cax * daz;
244        double cxdz = cax * day - cay * dax;
245
246        double da = dax*dax + day*day + daz*daz;
247        double ca = cax*cax + cay*cay + caz*caz;
248        double ba = bax*bax + bay*bay + baz*baz;
249
250        double nx = da * bxcx + ca * dxbx + ba * cxdx;
251        double ny = da * bxcy + ca * dxby + ba * cxdy;
252        double nz = da * bxcz + ca * dxbz + ba * cxdz;
253
254        double det = 2 * ( bax * cay * daz + bay * caz * dax + baz * cax * day
255                    - bax * caz * day - bay * cax * daz - baz * cay * dax);
256
257        sqrad = (nx*nx + ny*ny + nz*nz) / (det * det);
258        radius = sqrt(sqrad);
259        circumcentre->x = vertex[0]->x + nx / det;
260        circumcentre->y = vertex[0]->y + ny / det;
261        circumcentre->z = vertex[0]->z + nz / det;
262       
263        }
264        else
265        if(dimension == 2)
266        {
267/*
268        double els[9];
269        els[0] = vertex[1]->x - vertex[0]->x;
270        els[3] = vertex[1]->y - vertex[0]->y;
271        els[6] = vertex[1]->z - vertex[0]->z;
272        els[1] = vertex[2]->x - vertex[0]->x;
273        els[4] = vertex[2]->y - vertex[0]->y;
274        els[7] = vertex[2]->z - vertex[0]->z;
275        els[2] = els[3]*els[7] - els[6]*els[4];
276        els[5] = els[6]*els[1] - els[0]*els[7];
277        els[8] = els[0]*els[4] - els[3]*els[1];
278        Matrix M(3,3,els);
279
280        double vavg[3];
281        vavg[0] = (     els[0]*(vertex[1]->x + vertex[0]->x) +
282                                els[3]*(vertex[1]->y + vertex[0]->y) +
283                                els[6]*(vertex[1]->z + vertex[0]->z)
284                                )/2;
285        vavg[1] = (     els[1]*(vertex[2]->x + vertex[0]->x) +
286                                els[4]*(vertex[2]->y + vertex[0]->y) +
287                                els[7]*(vertex[2]->z + vertex[0]->z)
288                                )/2;
289        vavg[2] =       els[2]*vertex[0]->x +
290                                els[5]*vertex[0]->y +
291                                els[8]*vertex[0]->z;
292        Vector V(3,vavg);
293
294        M.invert();
295        Vector ResX = M*V;
296       
297        circumcentre->x = ResX(1);
298        circumcentre->y = ResX(2);
299        circumcentre->z = ResX(3);
300        double r1 = circumcentre->x - vertex[0]->x;
301        double r2 = circumcentre->y - vertex[0]->y;
302        double r3 = circumcentre->z - vertex[0]->z;
303        sqrad = r1*r1 + r2*r2 + r3*r3;
304        radius = sqrt(sqrad);
305*/     
306        }
307        //WARN 'dim=1' case missing
308}
309
310
311WsubMesh* Wsimplex::cavity(Vector *proxy)
312{
313        if(deleted)
314                return (WsubMesh*)0x00000001;
315        if(inCircle(proxy))
316        {
317                deleted = 1;
318
319                //WARN out of chain here, or undertaker later
320
321                //into deleted chain
322                WsubMesh* collector = new WsubMesh();
323               
324                collector->setInsider(this);
325
326                for(int i=0;i<dimension+1;i++)
327                {
328                        WsubMesh* a = NULL;
329                        if(neighbour[i])
330                        {
331                                a = neighbour[i]->cavity(proxy);
332                        }
333//else it is outer space, notice a == NULL then
334                        if(a)
335                        {
336//WARNING n[i] returned, therefore all his neighbours are dead,
337// he will never be called again, but cannot be deleted yet!!!
338// 0x1 return indicates a tetra earlier in call chain
339                                if(a != (WsubMesh*)0x00000001)
340                                {
341//                                      delete neighbour[i];
342                                        collector->merge(a);
343                                }
344                        }
345                        else
346                        {
347                                Wsimplex* n;
348                                if(dimension == 3)
349                                        n = new Wsimplex(
350                                        vertex[(i==0)?1:0],
351                                        vertex[(i<2)?2:1],
352                                        vertex[(i<3)?3:2]);
353                                else
354                                        n= new Wsimplex(
355                                        vertex[(i==0)?1:0],
356                                        vertex[(i==2)?1:2]);
357                                n->behind = neighbour[i];
358                                collector->add(n);
359                        }
360
361                }
362
363                return collector;
364               
365        }
366        else
367                return NULL;
368}
369
370void Wsimplex::helloNeighbour(Wsimplex* theguynext)
371{
372        Vector** ott = theguynext->getVertices();
373        if(!(vertex[0] == ott[0] || vertex[0] == ott[1]))
374        {
375                neighbour[0] = theguynext;
376                return;
377        }
378       
379        if(!(vertex[dimension] == ott[dimension] || vertex[dimension] == ott[dimension-1]))
380        {
381                neighbour[dimension] = theguynext;
382                return;
383        }
384
385        for(int t=1;t<dimension;t++)
386                if(!(vertex[t] == ott[t] || vertex[t] == ott[t+1] || vertex[t] == ott[t-1]))
387                {
388                        neighbour[t] = theguynext;
389                        return;
390                }
391        printf("ERROR: helloNeighbour found no matching nodes!\n");
392}
393
394int Wsimplex::getDimension()
395{
396        return dimension;
397}
398
399Vector** Wsimplex::getVertices()
400{
401        return vertex;
402}
403
404Wsimplex** Wsimplex::getNeighbours()
405{
406        return neighbour;
407}
408
409
410Wsimplex::Wsimplex(Vector *q, Vector *w, Vector *e, Vector *r)
411{
412
413//WARNING position is not checked for qualit
414//quality NOT SET!!
415//this constructor was marked obsolete
416//but it is used for constructing from file
417
418        dimension = 3;
419        deleted = 0;
420        for(int i=0;i<4;i++)
421                neighbour[i] = NULL;
422        vertex[0] = q;
423        vertex[1] = w;
424        vertex[2] = e;
425        vertex[3] = r;
426        circumcentre = new Vector(0,0,0);
427        updateCentre();
428}
429
430
431void Wsimplex::print()
432{
433//      cout << "Wsimplex:\n";
434//      for(int i=0;i<=dimension;i++)
435//              vertex[i]->print();
436}
437
438int Wsimplex::adjacentToNothing()
439{
440        int c=0;
441        for(int y=0;y<=dimension;y++)
442                if(!neighbour[y])
443                        c++;
444        return c;
445}
446
447int Wsimplex::verifyDelaunay(Vector* proxy)
448{
449        if(proxy)
450        {
451                if(inSmallSphere(proxy))
452                        return 0;
453                else
454                        return 1;
455        }
456        else
457        return 1;
458
459}
460
461int Wsimplex::inSmallSphere(Vector *guess)
462{
463        Vector xx = (*circumcentre) - (*guess);
464
465        return sqrad-0.1 > xx.norm2();
466}
467
468void Wsimplex::cutlinkback(Wbucket *container)
469{
470        if(!backlist.empty())
471        {
472                ::std::vector<Wbucket*>::iterator a = backlist.begin();
473                while(a != backlist.end())
474                {
475                        if(*a == container)
476                        {
477                                backlist.erase(a);
478                                break;
479                        }
480                        else
481                                a++;
482                }
483        }
484}
485
486Vector* Wsimplex::getInteriorPoint()
487{
488        double ix = 0, iy = 0, iz = 0;
489        for(int c=0;c<=dimension;c++)
490        {
491                ix += vertex[c]->x;
492                iy += vertex[c]->y;
493                iz += vertex[c]->z;
494        }
495        return new Vector(ix/(double)(dimension+1),iy/(double)(dimension+1),iz/(double)(dimension+1));
496}
497
498int Wsimplex::nearEdge()
499{
500//WARNING no id in Vector, should identify boundary nodes using geometry
501        return 0;
502
503}
504
505Wsimplex* Wsimplex::getNeighboursFront(int which)
506{
507        if(neighbour[which])
508                return neighbour[which]->getFront();
509        else
510                return 0x0;
511}
512
513Wsimplex* Wsimplex::getFront()
514{
515        return inFrontOf;
516}
517
518Wsimplex* Wsimplex::getBehind()
519{
520        return behind;
521}
522
523int Wsimplex::isComplete()
524{
525        Wsimplex** a = neighbour;
526        int com = 0;
527        for(int c=0; c<dimension; c++)
528                com = com && *(a++);
529        return com;
530}
531
532int Wsimplex::join(Wsimplex *other)
533{
534        Vector** ott = other->getVertices();
535
536        if(dimension == 1)
537        {
538                if(vertex[0] == ott[0])
539                {
540                        neighbour[1] = other;
541                        other->neighbour[1] = this;
542                }
543                else
544                        if(vertex[0] == ott[1])
545                        {
546                                neighbour[1] = other;
547                                other->neighbour[0] = this;
548                        }
549                else
550                        if(vertex[1] == ott[0])
551                        {
552                                neighbour[0] = other;
553                                other->neighbour[1] = this;
554                        }
555                else
556                        if(vertex[1] == ott[1])
557                        {
558                                neighbour[0] = other;
559                                other->neighbour[0] = this;
560                        }
561                return ((neighbour[0] && neighbour[1])?1:0) |
562                        ((other->neighbour[0] && other->neighbour[1])?2:0);
563        }
564
565        if(vertex[0] == ott[0])
566        {
567                if(vertex[1] == ott[1])
568                {
569                        neighbour[2] = other;
570                        other->neighbour[2] = this;
571                }
572                else
573                        if(vertex[1] == ott[2])
574                        {
575                                neighbour[2] = other;
576                                other->neighbour[1] = this;
577                        }
578                        else
579                                if(vertex[2] == ott[1])
580                                {
581                                        neighbour[1] = other;
582                                        other->neighbour[2] = this;
583                                }
584                                else
585                                        if(vertex[2] == ott[2])
586                                        {
587                                                neighbour[1] = other;
588                                                other->neighbour[1] = this;
589                                        }
590        }
591        else
592                if(vertex[0] == ott[1])
593                {
594                        if(vertex[1] == ott[2])
595                        {
596                                neighbour[2] = other;
597                                other->neighbour[0] = this;
598                        }
599                        else
600                                if(vertex[2] == ott[2])
601                                {
602                                        neighbour[1] = other;
603                                        other->neighbour[0] = this;
604                                }
605                }
606                else
607                        if(vertex[1] == ott[0])
608                        {
609                                if(vertex[2] == ott[1])
610                                {
611                                neighbour[0] = other;
612                                other->neighbour[2] = this;
613                                }
614                                else
615                                        if(vertex[2] == ott[2])
616                                        {
617                                        neighbour[0] = other;
618                                        other->neighbour[1] = this;
619                                        }
620                        }
621                        else
622                                if((vertex[1] == ott[1]) && (vertex[2] == ott[2]))
623                                        {
624                                        neighbour[0] = other;
625                                        other->neighbour[0] = this;
626                                        }
627                                else
628                                        return -1; //join impossible;
629//check if complete, return proper value
630        return ((neighbour[0] && neighbour[1] && neighbour[2])?1:0) |
631                        ((other->neighbour[0] && other->neighbour[1] && other->neighbour[2])?2:0);
632}
633
Note: See TracBrowser for help on using the repository browser.