← Back to branch summary

~ubuntu-branches/debian/sid/3depict/sid

« back to all changes in this revision

Viewing changes to src/isoSurface.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Sylvestre Ledru, D Haley, Sylvestre Ledru
  • Date: 2011-08-15 23:35:55 UTC
  • mfrom: (1.2.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20110815233555-5c0j9oklfkrc733q
Tags: 0.0.7-1
[ D Haley ]
* Update to new upstream version 
* Include manual in tarball
* Include upstream patches for correcting density profile computation
  and refresh branch computation

[ Sylvestre Ledru ]
* Switch to dpkg-source 3.0 (quilt) format

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/isoSurface.cpp
19
19
#include "isoSurface.h"
20
20
 
21
21
#include <limits>
 
22
#include <map>
 
23
 
 
24
 
 
25
#ifdef DEBUG
 
26
template<class T1, class T2>
 
27
bool mapUnique(const std::map<T1,T2> &m) 
 
28
{
 
29
        vector<T2> vec;
 
30
        vec.reserve(m.size());
 
31
 
 
32
        for(typename std::map<T1,T2>::const_iterator it = m.begin(); it!=m.end(); it++)
 
33
        {
 
34
                if(std::find(vec.begin(),vec.end(),it->second) != vec.end())
 
35
                        return false;
 
36
 
 
37
                vec.push_back(it->second);
 
38
 
 
39
        }
 
40
 
 
41
        return true;
 
42
}
 
43
#endif
 
44
 
 
45
 
 
46
//input vector "vec" must be sorted and unique
 
47
template<class T>
 
48
void removeElements( const std::vector<size_t> &elems,std::vector<T> &vec)
 
49
{
 
50
        if(vec.empty() || elems.empty())
 
51
                return;
 
52
 
 
53
        if(vec.size() == elems.size())
 
54
        {
 
55
                vec.clear();
 
56
                return;
 
57
        }
 
58
 
 
59
        size_t offset=vec.size()-1;     
 
60
        //Run backwards, swapping out 
 
61
        for(size_t ui=elems.size();ui--;)
 
62
        {
 
63
                ASSERT(ui <= offset);
 
64
                std::swap(vec[elems[ui]],vec[offset]);
 
65
                offset--;
 
66
        }
 
67
        vec.resize(offset+1);
 
68
}
 
69
 
 
70
 
 
71
 
 
72
 
22
73
void TriangleWithVertexNorm::computeACWNormal(Point3D &n) const
23
74
{
24
75
        Point3D a,b;
30
81
        n.normalise();
31
82
}
32
83
 
 
84
void TriangleWithVertexNorm::safeComputeACWNormal(Point3D &n) const
 
85
{
 
86
        Point3D a,b;
 
87
 
 
88
        a = p[0]-p[1];
 
89
        b = p[0]-p[2];
 
90
 
 
91
        n=a.crossProd(b);
 
92
        if(n.sqrMag() < sqrt(std::numeric_limits<float>::epsilon()) )
 
93
                n=Point3D(0,0,1);       
 
94
        else
 
95
                n.normalise();
 
96
 
 
97
 
 
98
 
 
99
}
 
100
 
 
101
float TriangleWithVertexNorm::computeArea() const
 
102
{
 
103
        Point3D a,b;
 
104
 
 
105
        a = p[0]-p[1];
 
106
        b = p[0]-p[2];
 
107
 
 
108
        return  a.crossProd(b).sqrMag();
 
109
 
 
110
}
 
111
 
33
112
bool TriangleWithVertexNorm::isDegenerate() const
34
113
{
35
114
        return (p[0].sqrDist(p[1]) < std::numeric_limits<float>::epsilon() ||
37
116
        p[2].sqrDist(p[1]) < std::numeric_limits<float>::epsilon());
38
117
}
39
118
 
 
119
void TriangleWithVertexNorm::getCentroid(Point3D &c) const
 
120
{
 
121
        c=p[0];
 
122
        c+=p[1];
 
123
        c+=p[2];
 
124
 
 
125
        c*=1.0f/3.0f;
 
126
}
 
127
 
40
128
//This code is a modified version of the following:
41
129
//==============
42
130
// Marching Cubes Example Program 
354
442
        {0.0, 0.0, 1.0},{1.0, 0.0, 1.0},{1.0, 1.0, 1.0},{0.0, 1.0, 1.0}
355
443
};
356
444
 
 
445
//The deltas for the vertex offsets.  This is the int version of a2fVertexOffset
 
446
const unsigned int VERTEX_OFFSET[8][3] = 
 
447
{
 
448
        {0, 0, 0},{1, 0, 0},{1, 1, 0},{0, 1, 0},
 
449
        {0, 0, 1},{1, 0, 1},{1, 1, 1},{0, 1, 1}
 
450
};
 
451
 
357
452
//a2iEdgeConnection lists the index of the endpoint vertices for each of the 12 edges of the cube
358
453
static const int a2iEdgeConnection[12][2] = 
359
454
{
370
465
        {0.0, 0.0, 1.0},{0.0, 0.0, 1.0},{ 0.0, 0.0, 1.0},{0.0,  0.0, 1.0}
371
466
};
372
467
 
 
468
//Mapping between edges as defined in a2iEdgeConenction, (and thus implicitly in edge table)
 
469
//and the voxels.h definition
 
470
int edgeRemap[12]  ={ 0,6,1,4,
 
471
                      2,7,3,5,
 
472
                      8,10,11,9};
 
473
 
 
474
 
 
475
 
373
476
 
374
477
//fGetOffset finds the approximate point of intersection of the surface
375
478
// between two points with the values fValue1 and fValue2
383
486
        return (fValueDesired - fValue1)/fDelta;
384
487
}
385
488
 
 
489
 
386
490
//vMarchingCubes iterates over the entire dataset, calling vMarchCube on each cube
387
491
void marchingCubes(const Voxels<float> &v,float isoValue, vector<TriangleWithVertexNorm> &tVec)
388
492
{
390
494
        v.getSize(nx,ny,nz);
391
495
 
392
496
        ASSERT(nx > 1 && ny>1 && nz>1);
 
497
 
 
498
        //Don't try to isosurface a any volume with a unitary dimension.
 
499
        if(nx ==1 || ny ==1 || nz == 1)
 
500
                return;
 
501
 
393
502
        Point3D gridSpacing;
394
503
        gridSpacing=v.getPitch();
395
504
 
397
506
        BoundCube boundC;
398
507
        boundC.setBounds(v.getMinBounds(),v.getMaxBounds());
399
508
#endif
400
 
        float vertexVal[8];
 
509
 
 
510
        vector<TriangleWithIndexedVertices> indexedTriVec;
 
511
 
401
512
        //Loop over the vertexs, with the mesh such that the 
402
513
        //nominally cube centres are now on a grid that is dual
403
514
        //to the original grid (excluding the external boundary of course)
 
515
#pragma omp parallel for
404
516
        for(size_t iX = 0; iX < nx-1; iX++)
405
 
        for(size_t iY = 0; iY < ny-1; iY++)
406
 
        for(size_t iZ = 0; iZ < nz-1; iZ++)
407
 
        {
408
 
                int iEdgeFlags,iFlagIndex;
409
 
                iEdgeFlags=iFlagIndex=0;
410
 
                Point3D position;
411
 
                //Lower left corner of cell for dual grid
412
 
                position=v.getPoint(iX,iY,iZ) + gridSpacing*0.5;
413
 
 
414
 
                //This must match a2fVertexOffset
415
 
                vertexVal[0] = v.getData(iX,iY,iZ);
416
 
                vertexVal[1] = v.getData(iX+1,iY,iZ);
417
 
                vertexVal[2] = v.getData(iX+1,iY+1,iZ);
418
 
                vertexVal[3] = v.getData(iX,iY+1,iZ);
419
 
 
420
 
                vertexVal[4] = v.getData(iX,iY,iZ+1);
421
 
                vertexVal[5] = v.getData(iX+1,iY,iZ+1);
422
 
                vertexVal[6] = v.getData(iX+1,iY+1,iZ+1);
423
 
                vertexVal[7] = v.getData(iX,iY+1,iZ+1);
424
 
 
425
 
                //Find which vertices are inside of the surface and which are outside
426
 
                for(int iVertexTest = 0; iVertexTest < 8; iVertexTest++)
427
 
                {
428
 
                        if(vertexVal[iVertexTest] <= isoValue) 
429
 
                                iFlagIndex |= 1<<iVertexTest;
430
 
                }
431
 
                //Find which edges are intersected by the surface
432
 
                iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
433
 
                if(!iEdgeFlags) 
 
517
        {
 
518
                int iEdgeFlags,iFlagIndex;      
 
519
                for(size_t iY = 0; iY < ny-1; iY++)
 
520
                {
 
521
                for(size_t iZ = 0; iZ < nz-1; iZ++)
 
522
                {
 
523
                        iEdgeFlags=iFlagIndex=0;
 
524
                        Point3D position;
 
525
                        //Lower left corner of cell for dual grid
 
526
                        position=v.getPoint(iX,iY,iZ) + gridSpacing*0.5;
 
527
 
 
528
                        //Find which vertices are inside of the surface and which are outside
 
529
                        for(int iVertexTest = 0; iVertexTest < 8; iVertexTest++)
 
530
                        {
 
531
                                float f;
 
532
                                f=v.getData(iX+VERTEX_OFFSET[iVertexTest][0],
 
533
                                                iY+VERTEX_OFFSET[iVertexTest][1],
 
534
                                                iZ+VERTEX_OFFSET[iVertexTest][2]);
 
535
 
 
536
                                //Compute posiion in triangle and edge connection
 
537
                                //tables
 
538
                                if(f <= isoValue) 
 
539
                                        iFlagIndex |= 1<<iVertexTest;
 
540
                        }
 
541
                        //Find which edges are intersected by the surface
 
542
                        iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
 
543
                        if(!iEdgeFlags) 
 
544
                                continue;
 
545
 
 
546
 
 
547
                        //Find the point of intersection of the surface with each edge
 
548
                        //Then find the normal to the surface at those points
 
549
                        size_t asEdgeVertex[12];
 
550
                        for(int iEdge = 0; iEdge < 12; iEdge++)
 
551
                        {
 
552
                                //if there is an intersection on this edge
 
553
                                if(iEdgeFlags & (1<<iEdge))
 
554
                                {
 
555
                                        //Store a  reference to the vertex position
 
556
                                        asEdgeVertex[iEdge] = v.getEdgeIndex(iX,iY,iZ,edgeRemap[iEdge]);
 
557
                                }
 
558
                        }
 
559
 
 
560
                
 
561
                        //Store the triangles that were found.  There can be up to five per cube; 
 
562
                        //these are listed as triplets in the connection table
 
563
                        for(int iTriangle = 0; iTriangle < 5; iTriangle++)
 
564
                        {
 
565
                                //Check to see if this triplet is a valid triangle
 
566
                                if(a2iTriangleConnectionTable[iFlagIndex][3*iTriangle] < 0)
 
567
                                        break;
 
568
 
 
569
                                TriangleWithIndexedVertices t;
 
570
                                for(int iCorner = 0; iCorner < 3; iCorner++)
 
571
                                {
 
572
                                        int iVertex;
 
573
                                        iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
 
574
                                        //we should only be accessing an edge if the edge was set.
 
575
                                        ASSERT((1 << iVertex) & iEdgeFlags);
 
576
 
 
577
 
 
578
 
 
579
 
 
580
                                        t.p[iCorner] = asEdgeVertex[iVertex];
 
581
                                }
 
582
#pragma omp critical
 
583
                                indexedTriVec.push_back(t);
 
584
                        }
 
585
                }
 
586
                }
 
587
        }
 
588
 
 
589
        //---------
 
590
 
 
591
        //OK, so we set up a network of triangle vertices which consist of edge numbers
 
592
        //and a list of triangles which are the triplets of those numbers.
 
593
        //Convert this into a 
 
594
 
 
595
        //Map that contains a list of triangles in the vector  "indexedTriVec"
 
596
        // which share a particular edge (or rather, position along that edge)
 
597
        // First element is the edge index. Second element is the list of triangles
 
598
        // who share this edge.
 
599
        //
 
600
        // We will need this when we do the vertex normals.
 
601
        std::map<size_t,list<size_t> > edgeTriMap;
 
602
#pragma omp parallel for 
 
603
        for(size_t ui=0;ui<indexedTriVec.size(); ui++)  
 
604
        {
 
605
                for(size_t uj=0;uj<3;uj++)
 
606
                {
 
607
                        map<size_t,list<size_t> >::iterator it;
 
608
                        it = edgeTriMap.find((indexedTriVec[ui].p[uj] >>2 ) << 2);
 
609
                        if(it == edgeTriMap.end())
 
610
                        {
 
611
 
 
612
                                list<size_t> seedList;
 
613
                                seedList.push_back(ui);
 
614
#pragma omp critical
 
615
                                edgeTriMap.insert(
 
616
                                        make_pair((indexedTriVec[ui].p[uj] >>2 ) << 2,seedList));
 
617
                                                
 
618
                        }
 
619
                        else
 
620
                        {
 
621
                                it->second.push_back(ui);
 
622
                        }
 
623
                }
 
624
 
 
625
        }
 
626
 
 
627
 
 
628
        //Generate the position points for each edge
 
629
        map<size_t,Point3D> pointMap;
 
630
        for(map<size_t,list<size_t> >::iterator it=edgeTriMap.begin();
 
631
                it!=edgeTriMap.end(); ++it)
 
632
        {
 
633
                Point3D low,high,voxelFrameIntersection; 
 
634
                float lowF,highF,alpha;
 
635
 
 
636
                if(pointMap.find((it->first>>2)<<2) != pointMap.end())
434
637
                        continue;
435
 
                //Find the point of intersection of the surface with each edge
436
 
                //Then find the normal to the surface at those points
437
 
                Point3D asEdgeVertex[12];
438
 
                for(int iEdge = 0; iEdge < 12; iEdge++)
439
 
                {
440
 
                        //if there is an intersection on this edge
441
 
                        if(iEdgeFlags & (1<<iEdge))
442
 
                        {
443
 
                                float fOffset;
444
 
                                fOffset =0.5;// interpLinear(vertexVal[ a2iEdgeConnection[iEdge][0] ], 
445
 
                                        //      vertexVal[ a2iEdgeConnection[iEdge][1] ], isoValue);
446
 
 
447
 
                                ASSERT(fOffset >= -0.01 && fOffset < 1.01 );
448
 
 
449
 
                                asEdgeVertex[iEdge][0]= position[0] + 
450
 
                                        (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][0]  +  
451
 
                                                fOffset * a2fEdgeDirection[iEdge][0]) * gridSpacing[0];
452
 
 
453
 
                                asEdgeVertex[iEdge][1] = position[1]+ 
454
 
                                        (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][1]  +  
455
 
                                         fOffset * a2fEdgeDirection[iEdge][1]) * gridSpacing[1];
456
 
                                asEdgeVertex[iEdge][2] = position[2] + 
457
 
                                        (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][2]  +  
458
 
                                         fOffset * a2fEdgeDirection[iEdge][2]) * gridSpacing[2];
459
 
 
460
 
 
461
 
                                ASSERT(boundC.containsPt(asEdgeVertex[iEdge] ));
462
 
                                //We will have to set normal later
463
 
                        }
464
 
                }
465
 
 
466
 
        
467
 
                //Store the triangles that were found.  There can be up to five per cube
468
 
                for(int iTriangle = 0; iTriangle < 5; iTriangle++)
469
 
                {
470
 
                        if(a2iTriangleConnectionTable[iFlagIndex][3*iTriangle] < 0)
471
 
                                break;
472
 
 
473
 
                        TriangleWithVertexNorm t;
474
 
                        for(int iCorner = 0; iCorner < 3; iCorner++)
475
 
                        {
476
 
                                int iVertex;
477
 
                                iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
478
 
 
479
 
                                //we should only be accessing an edge if the edge was set.
480
 
                                ASSERT((1 << iVertex) & iEdgeFlags);
481
 
                                t.p[iCorner] = asEdgeVertex[iVertex];
482
 
                                //The mesh must be within the grid
483
 
                                ASSERT(boundC.containsPt(t.p[iCorner]));
484
 
                        }
485
 
 
486
 
                        if(!t.isDegenerate())
487
 
                        {
488
 
                                Point3D p;
489
 
                                t.computeACWNormal(p);
490
 
                                
491
 
                                //Assign pure normals,
492
 
                                for(unsigned int ui=0;ui<3;ui++)
493
 
                                        t.normal[ui]=p;
494
 
 
495
 
                                tVec.push_back(t);
496
 
                        }
497
 
                }
498
 
        }
499
 
 
500
 
 
501
 
        
 
638
 
 
639
                //Low/high sides of edge's scalar values
 
640
                v.getEdgeEndApproxVals(it->first,lowF,highF);
 
641
 
 
642
 
 
643
                //Get the edge's low and high end node positions
 
644
                v.getEdgeEnds((it->first>>2)<<2,low,high);
 
645
                
 
646
                //OK, now we have that, lets use the values to "lever" the 
 
647
                //solution point note node locations for isosurface 
 
648
                if(fabs(highF-lowF) < sqrt(std::numeric_limits<float>::epsilon()))
 
649
                {
 
650
                        //Prevent divide by zero
 
651
                        voxelFrameIntersection=(low+high)*0.5;
 
652
                }
 
653
                else
 
654
                {
 
655
                        //interpolate
 
656
                        alpha= (isoValue- lowF) / (highF- lowF);
 
657
                        voxelFrameIntersection=low + (high-low)*alpha;
 
658
                }
 
659
 
 
660
                        
 
661
                pointMap.insert(make_pair((it->first>>2)<<2,voxelFrameIntersection));
 
662
        }       
 
663
 
 
664
 
 
665
        tVec.resize(indexedTriVec.size());
 
666
        vector<size_t> popTris;
 
667
        //Set all triangle verticies
 
668
        #pragma omp parallel for
 
669
        for(size_t ui=0;ui<tVec.size();ui++)
 
670
        {
 
671
        
 
672
                //Do a lookup of the edge point locations
 
673
                // by mapping the edge indicies to the edge poitns 
 
674
                for(int uj=0;uj<3;uj++)
 
675
                        tVec[ui].p[uj] = pointMap.at((indexedTriVec[ui].p[uj]>>2)<<2);
 
676
 
 
677
                if(tVec[ui].isDegenerate())
 
678
                {
 
679
#pragma omp critical
 
680
                        popTris.push_back(ui);
 
681
                }
 
682
        }
 
683
        
 
684
 
 
685
        //Remove any degenerate triangles from both indexed and real maps
 
686
        removeElements(popTris,tVec);
 
687
        removeElements(popTris,indexedTriVec);
 
688
 
 
689
        if(!tVec.size())
 
690
                return;
 
691
        
 
692
        //set all triangle edge normals by inverse face area weighting.
 
693
        // The idea is that big triangles don't affect the normal at the pooint
 
694
        // as they are quite delocalised. Little triangles affect it mroe.
 
695
        // This is entirely empirical
 
696
        // ----
 
697
        
 
698
        vector<Point3D> origNormal;
 
699
        origNormal.resize(indexedTriVec.size());
 
700
        #pragma omp parallel for
 
701
        for(size_t ui=0;ui<indexedTriVec.size();ui++)
 
702
                tVec[ui].safeComputeACWNormal(origNormal[ui]);
 
703
 
 
704
 
 
705
        //Loop over all vertices again, then asssign a weighted
 
706
        //normal in place of the original normal. Lets cheat and
 
707
        //re-use "pointmap".
 
708
        // --
 
709
        for(map<size_t,Point3D>::iterator it=pointMap.begin(); 
 
710
                                        it!=pointMap.end();++it)
 
711
                it->second=Point3D(0,0,0);
 
712
 
 
713
        //Construct the shared normals
 
714
        float smallNum=sqrt(std::numeric_limits<float>::epsilon());
 
715
        for(size_t ui=0;ui<indexedTriVec.size();ui++)
 
716
        {
 
717
                //For eachvertex in our current triangle
 
718
                //update the weight mapping
 
719
                float weight;
 
720
                weight=tVec[ui].computeArea();
 
721
 
 
722
                if(weight > smallNum)
 
723
                {
 
724
                        for(int uj=0;uj<3;uj++)
 
725
                                pointMap.at((indexedTriVec[ui].p[uj]>>2) <<2)+=origNormal[ui]*weight;
 
726
                }
 
727
        }
 
728
 
 
729
 
 
730
        //re-normalise normals
 
731
        for(map<size_t,Point3D>::iterator it=pointMap.begin(); 
 
732
                                        it!=pointMap.end();++it)
 
733
        {
 
734
                if(it->second.sqrMag() > smallNum)
 
735
                        it->second.normalise();
 
736
                else
 
737
                        it->second=Point3D(0,0,1);
 
738
        }
 
739
 
 
740
        //assign these normals to the vertices of each triangle
 
741
        #pragma omp parallel for
 
742
        for(size_t ui=0;ui<indexedTriVec.size();ui++)
 
743
        {
 
744
                for(unsigned int uj=0;uj<3;uj++)
 
745
                        tVec[ui].normal[uj] = pointMap.at((indexedTriVec[ui].p[uj]>>2)<<2);
 
746
        }       
 
747
        // --
 
748
        
 
749
        // ----
 
750
 
 
751
 
 
752
        //TODO: We could use something like triStripper 
 
753
        // to optimise this, rather than just build the raw triangles
 
754
        // http://users.telenet.be/tfautre/softdev/tristripper/
 
755
 
502
756
}