GeoTessGrid.h

Go to the documentation of this file.

//- ****************************************************************************
//- 
//- Copyright 2009 Sandia Corporation. Under the terms of Contract
//- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
//- retains certain rights in this software.
//- 
//- BSD Open Source License.
//- All rights reserved.
//- 
//- Redistribution and use in source and binary forms, with or without
//- modification, are permitted provided that the following conditions are met:
//- 
//-    * Redistributions of source code must retain the above copyright notice,
//-      this list of conditions and the following disclaimer.
//-    * Redistributions in binary form must reproduce the above copyright
//-      notice, this list of conditions and the following disclaimer in the
//-      documentation and/or other materials provided with the distribution.
//-    * Neither the name of Sandia National Laboratories nor the names of its
//-      contributors may be used to endorse or promote products derived from
//-      this software without specific prior written permission.
//- 
//- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
//- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
//- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
//- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
//- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
//- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
//- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
//- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
//- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
//- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//- POSSIBILITY OF SUCH DAMAGE.
//-
//- ****************************************************************************
 
#ifndef GEOTESSGRID_OBJECT_H
#define GEOTESSGRID_OBJECT_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
 
#include <iostream>
#include <string>
#include <fstream>
#include <vector>
#include <set>
#include <list>
#include <set>
#include <sstream>
 
// use standard library objects
using namespace std;
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "CPPUtils.h"
#include "GeoTessUtils.h"
#include "ArrayReuse.h"
#include "GeoTessException.h"
#include "GeoTessOptimizationType.h"
 
// **** _BEGIN GEOTESS NAMESPACE_ **********************************************
 
namespace geotess
{
 
// **** _FORWARD REFERENCES_ ***************************************************
 
class IFStreamAscii;
class IFStreamBinary;
 
// **** _CLASS DEFINITION_ *****************************************************
 
struct Edge
{
    int vj;
 
    int vk;
 
    int cornerj;
 
    int tLeft;
 
    int tRight;
 
    double normal[3]; //
 
    Edge* next; //
};
 
class GEOTESS_EXP_IMP GeoTessGrid
{
protected:
 
    double** vertices;
 
    int nVertices;
 
    int** triangles;
 
    int nTriangles;
 
    int** levels;
 
    int nLevels;
 
    int** tessellations;
 
    int nTessellations;
 
    void initialize();
 
    GeoTessGrid(const string& gid) :
        vertices(NULL), nVertices(0),
        triangles(NULL), nTriangles(0),
        levels(NULL), nLevels(0),
        tessellations(NULL), nTessellations(0),
        descendants(NULL), gridID(gid),
        gridInputFile("null"), gridOutputFile("null"), gridSoftwareVersion(""),
        gridGenerationDate(""), refCount(0)
    {
    }
 
private:
 
    int* descendants;
 
    vector<vector<Edge*> > edgeList;
 
    mutable vector<vector<Edge*> > spokeList;
 
    mutable vector<double*> circumCenters;
 
    string gridID;
 
    string gridInputFile;
 
    string gridOutputFile;
 
    string gridSoftwareVersion;
 
    string gridGenerationDate;
 
    int refCount;
 
    vector< vector< vector<int> > > vtxTriangles;
 
    mutable vector< set<int> > connectedVertices;
 
    GeoTessGrid* loadGridAscii(const string& inputFile);
 
    GeoTessGrid* loadGridAscii(IFStreamAscii& input);
 
    GeoTessGrid* loadGridBinary(const string& inputFile);
 
    virtual GeoTessGrid* loadGridBinary(IFStreamBinary& ifs);
 
    static void loadGridAsciiFront(IFStreamAscii& input, int& gridFileFormat,
            string& gridSWVersion, string& fileCreationDate,
            string& gridid, const string& grdInptFile);
 
    static void loadGridBinaryFront(IFStreamBinary& ifs, int& gridFileFormat,
            string& gridSWVersion, string& fileCreationDate,
            string& gridid, const string& grdInptFile);
 
    const set<int>& getVertexIndices(const int& level);
 
    void computeSpokeLists(const int& level) const;
 
public:
 
    GeoTessGrid() :
        vertices(NULL), nVertices(0),
        triangles(NULL), nTriangles(0),
        levels(NULL), nLevels(0),
        tessellations(NULL), nTessellations(0),
        descendants(NULL), gridID(""),
        gridInputFile("null"), gridOutputFile("null"), gridSoftwareVersion(""),
        gridGenerationDate(""), refCount(0)
    { }
 
    GeoTessGrid(IFStreamAscii& input) :
        vertices(NULL), triangles(NULL), levels(NULL), tessellations(NULL),
        descendants(NULL), gridID(""),
        gridInputFile("null"), gridOutputFile("null"), gridSoftwareVersion(""),
        gridGenerationDate(""), refCount(0)
    { loadGridAscii(input); }
 
    GeoTessGrid(IFStreamBinary& input) :
        vertices(NULL), triangles(NULL), levels(NULL), tessellations(NULL),
        descendants(NULL), gridID(""),
        gridInputFile("null"), gridOutputFile("null"), gridSoftwareVersion(""),
        gridGenerationDate(""), refCount(0)
    { loadGridBinary(input); }
 
    GeoTessGrid(GeoTessGrid &other);
 
    GeoTessGrid& operator= (const GeoTessGrid& other);
 
    LONG_INT getMemory()
    {
        LONG_INT memory = (LONG_INT)sizeof(GeoTessGrid);
 
        // double** vertices
        memory += nVertices * (LONG_INT)sizeof(double*) + nVertices * 3 * (LONG_INT)sizeof(double);
 
        // int** triangles
        memory += nTriangles * (LONG_INT)sizeof(int*) + nTriangles * 3 * (LONG_INT)sizeof(int);
 
        // int** levels
        memory += nLevels * (LONG_INT)sizeof(int*) + nLevels * 2 * (LONG_INT)sizeof(int);
 
        // int** tessellations
        memory += nTessellations * (LONG_INT)sizeof(int*) + nTessellations * 2 * (LONG_INT)sizeof(int);
 
        // int* descendants
        memory += nTriangles * (LONG_INT)sizeof(int);
 
        // vector<vector<Edge*> > edgeList : nTriangles x 3 array of Edges
        memory += (LONG_INT) (edgeList.capacity() * sizeof(vector<Edge*>)
                + edgeList.size() * 3 * (sizeof(Edge*) + sizeof(Edge)));
 
        // vector<vector<Edge*> > spokeList;
        memory += (LONG_INT) (spokeList.capacity() * sizeof(vector<Edge*>)
                + spokeList.size() * 3 * sizeof(Edge*));
 
        // vector<double*> circumCenters is nTriangles x 3 array of doubles
        memory += (LONG_INT) (circumCenters.capacity() * sizeof(double*));
        for (int i=0; i<(int)circumCenters.size(); ++i)
            if (circumCenters[i]) memory += 3 * (LONG_INT)sizeof(double);
 
        // An nLevels x nVertices x n array that stores the indices of the n
        // triangles of which each vertex is a member.
        // vector< vector< vector<int> > > vtxTriangles;
        memory += (LONG_INT) (vtxTriangles.capacity() * sizeof(vector< vector<int> >));
        for (int i=0; i<(int)vtxTriangles.size(); ++i)
        {
            memory += (LONG_INT) (vtxTriangles[i].capacity() * sizeof(vector<int>));
            for (int j=0; j<(int)vtxTriangles[i].size(); ++j)
                memory += (LONG_INT) (vtxTriangles[i][j].capacity() * sizeof(int));
        }
 
        // An nLevels x n array that stores the indexes of the n vertices
        // that are connected together by triangles on the corresponding
        // level.  Lazy evaluation is used to construct this array of sets.
        // vector< set<int> > connectedVertices;
        memory += (LONG_INT) (connectedVertices.capacity() * sizeof(set<int>));
        for (int i=0; i<(int)connectedVertices.size(); ++i)
            memory += (LONG_INT) (connectedVertices[i].size() * sizeof(int));
 
        // add memory requirements for all the string variables.
        memory += (LONG_INT) (gridID.length() + gridInputFile.length() + gridOutputFile.length()
                + gridSoftwareVersion.length() + gridGenerationDate.length());
 
        return memory;
    }
 
    GeoTessGrid* loadGrid(const string& inputFile);
 
    static string getGridID(const string& fileName);
 
    static bool isGeoTessGrid(const string& inputFile);
 
    bool                    operator == (const GeoTessGrid& g) const
                                                                            {   return (gridID == g.gridID); }
 
    bool                    operator != (const GeoTessGrid& g) const
                                                                            {   return !(*this == g); }
 
    const string& getGridID() const
    {
        return gridID;
    }
 
    void setGridSoftwareVersion(const string& swVersion) {gridSoftwareVersion = swVersion; };
 
    const string& getGridSoftwareVersion() const {return gridSoftwareVersion; };
 
    void setGridGenerationDate(const string& gridDate) {gridGenerationDate = gridDate; };
 
    const string& getGridGenerationDate() const {return gridGenerationDate; };
 
    void setGridInputFile(const string& gridFile) { gridInputFile = gridFile; }
 
    const string& getGridInputFile() const {return gridInputFile; }
 
    const string& getGridOutputFile() const { return gridOutputFile; }
 
    const double* getVertex(int vertex) const
    {
        return vertices[vertex];
    }
 
    int findClosestVertex(double* unit_vector, int tessId)
    {
        int* t = triangles[getTriangle(getFirstTriangle(tessId, 0), unit_vector)];
 
        int index = 0;
        double dot = GeoTessUtils::dot(unit_vector, vertices[t[0]]);
 
        double doti = GeoTessUtils::dot(unit_vector, vertices[t[1]]);
        if (doti > dot)
        {
            index = 1;
            dot = doti;
        }
 
        doti = GeoTessUtils::dot(unit_vector, vertices[t[2]]);
        if (doti > dot)
        {
            index = 2;
            dot = doti;
        }
 
        return t[index];
    }
 
    int getVertexIndex(int triangle, int corner) const
    { return triangles[triangle][corner]; }
 
    int getVertexIndex(int tessId, int level, int triangle, int corner) const
    {
        return triangles[levels[tessellations[tessId][0] + level][0] + triangle][corner];
    }
 
    int getVertexIndex(const double* u)
    {
        for (int i=nTessellations-1; i>=0; --i)
        {
            int vid = getVertexIndex(u, i);
            if (vid >= 0)
                return vid;
        }
        return -1;
    }
 
    int getVertexIndex(const double* u, int tessId)
    {
        int* t = triangles[getTriangle(getFirstTriangle(tessId, 0), u)];
 
        if (GeoTessUtils::dot(u, vertices[t[0]]) > cos(1e-7))
            return t[0];
        if (GeoTessUtils::dot(u, vertices[t[1]]) > cos(1e-7))
            return t[1];
        if (GeoTessUtils::dot(u, vertices[t[2]]) > cos(1e-7))
            return t[2];
        return -1;
    }
 
    double* getVertex(int tessId, int level, int triangle, int corner)
    {
        return vertices[triangles[levels[tessellations[tessId][0] + level][0] + triangle][corner]];
    }
 
    double const* const* getVertices() const { return vertices; }
 
    void getVerticesTopLevel(const int& tessellation, set<const double*>& vectors)
    {
        getVertices(tessellation, getNLevels(tessellation)-1, vectors);
    }
 
    void getVertices(const int& tessellation, const int& level, set<const double*>& vectors);
 
    const set<int>& getVertexIndices(const int& tessId, const int& level)
                                    { return getVertexIndices(getLevel(tessId, level)); }
 
    const set<int>& getVertexIndicesTopLevel(const int& tessId)
                                    { return getVertexIndices(tessellations[tessId][1] - 1); }
 
    int getNVertices() const { return nVertices; }
 
    int getNTessellations() const { return nTessellations; }
 
    int getNLevels(int tessellation) const
    { return tessellations[tessellation][1] - tessellations[tessellation][0]; }
 
    int getLevel(int tessellation, int i) const
    { return tessellations[tessellation][0] + i; }
 
    int getLastLevel(int tessellation) const
    { return tessellations[tessellation][1] - 1; }
 
    int getTopLevel(int tessellation) const
    { return tessellations[tessellation][1] - tessellations[tessellation][0] - 1; }
 
    int getNTriangles(int tessellation, int level) const
    {
        return levels[tessellations[tessellation][0] + level][1]
                                                              - levels[tessellations[tessellation][0] + level][0];
    }
 
    int getTriangle(int tessellation, int level, int i) const
    { return levels[tessellations[tessellation][0] + level][0] + i; }
 
    int getFirstTriangle(int tessellation, int level) const
    { return levels[tessellations[tessellation][0] + level][0]; }
 
    int getLastTriangle(int tessellation, int level) const
    { return levels[tessellations[tessellation][0] + level][1]-1; }
 
    int const* const* getTriangles() const { return triangles; }
 
    const int* getTriangleVertexIndexes(int triangleIndex) const
    { return triangles[triangleIndex]; }
 
    int getTriangleVertexIndex(int triangleIndex, int cornerIndex) const
    { return triangles[triangleIndex][cornerIndex]; }
 
    const double* getTriangleVertex(int triangleIndex, int cornerIndex) const
    { return vertices[triangles[triangleIndex][cornerIndex]]; }
 
 
    void getTriangleVertices(int triangle, double** triVrt)
    {
        int* corners = triangles[triangle];
        triVrt[0][0] = vertices[corners[0]][0];
        triVrt[0][1] = vertices[corners[0]][1];
        triVrt[0][2] = vertices[corners[0]][2];
        triVrt[1][0] = vertices[corners[1]][0];
        triVrt[1][1] = vertices[corners[1]][1];
        triVrt[1][2] = vertices[corners[1]][2];
        triVrt[2][0] = vertices[corners[2]][0];
        triVrt[2][1] = vertices[corners[2]][1];
        triVrt[2][2] = vertices[corners[2]][2];
    }
 
    const void computeCircumCenters()
    {
        if (circumCenters.empty()) circumCenters.resize(nTriangles, NULL);
 
        for (int triangle=0; triangle<nTriangles; ++triangle)
            if (circumCenters[triangle] == NULL)
            {
                int* corners = triangles[triangle];
                circumCenters[triangle] = GeoTessUtils::circumCenterPlus(vertices[corners[0]],
                        vertices[corners[1]], vertices[corners[2]]);
            }
    }
 
    const void computeCircumCenters(const int& level)
    {
        if (circumCenters.empty()) circumCenters.resize(nTriangles);
 
        if (circumCenters[levels[level][0]] == NULL)
        {
            for (int triangle=levels[level][0]; triangle<levels[level][1]; ++triangle)
            {
                int* corners = triangles[triangle];
                circumCenters[triangle] =  GeoTessUtils::circumCenterPlus(vertices[corners[0]],
                        vertices[corners[1]], vertices[corners[2]]);
            }
        }
    }
 
    const double* getCircumCenter(const int& triangle) const { return circumCenters[triangle]; }
 
    void getCircumCenter(const int& triangle, double* cc) const
    { cc[0]=circumCenters[triangle][0]; cc[1]=circumCenters[triangle][1]; cc[2]=circumCenters[triangle][2]; }
 
    int getNeighbor(const int& triangleIndex, const int& neighborIndex) const
    { return neighborIndex == 3 ? descendants[triangleIndex] : edgeList[triangleIndex][neighborIndex]->tLeft; }
 
    int getNeighbor(const int& tessellation, const int& level, const int& triangle, const int& side)
    { return getNeighbor(getTriangle(tessellation, level, triangle), side); }
 
    void getNeighbors(int triangleIndex, vector<int>& neighbors)
    {
        neighbors.clear();
        vector<Edge*> edges = edgeList[triangleIndex];
        neighbors.push_back(edges[0]->tLeft);
        neighbors.push_back(edges[1]->tLeft);
        neighbors.push_back(edges[2]->tLeft);
        neighbors.push_back(descendants[triangleIndex]);
    }
 
    void getNeighbors(const int& tessellation, const int& level, const int& triangle, vector<int>& neighbors)
    { getNeighbors(getTriangle(tessellation, level, triangle), neighbors); }
 
    int getNeighborIndex(const int& tid, const int& nid)
    {
        if (edgeList[tid][0]->tLeft == nid)
            return 0;
        if (edgeList[tid][1]->tLeft == nid)
            return 1;
        if (edgeList[tid][2]->tLeft == nid)
            return 2;
        return -1;
    }
 
    int const* getDescendants() const { return descendants; }
 
    int getDescendant(const int& triangle) const { return descendants[triangle]; }
 
    int getDescendant(const int& tessId, const int& level, const int& triangle) const
    { return descendants[getTriangle(tessId, level, triangle)]; }
 
    string toString();
 
    int getTriangle(int triangleIndex, const double* vector);
 
    const vector<vector<int> >& getVertexTriangles(const int& level) const { return vtxTriangles[level]; }
 
    const vector<int>& getVertexTriangles(const int& tessId, const int& level, const int& vertex) const
                            { return vtxTriangles[tessellations[tessId][0] + level][vertex]; }
 
    const vector<int>& getVertexTriangles(int tessId, int vertex) const
                            { return getVertexTriangles(tessId, getNLevels(tessId)-1, vertex); }
 
    //void clearVertexTriangles();
 
    void getVertexNeighborsOrdered(const int& tessId, const int& level, const int& vertex, vector<int>& v);
 
    void getVertexNeighbors(const int& tessId, const int& level, const int& vertex, set<int>& nbrs)
    { getVertexNeighbors(tessId, level, vertex, 1, nbrs); }
 
    void getVertexNeighbors(const int& tessId, const int& level, const int& vertex,
            const int& order, set<int>& nbrs);
 
    // All methods below this point are public but are not documented in the doxygen documentation.
    // These are methods that typical applications will never need to call.  They have to be
    // public because other classes in the GeoTess namespace need to access them.
    //
 
    GeoTessGrid(const GeoTessOptimizationType* opttype) :
        vertices(NULL), nVertices(0),
        triangles(NULL), nTriangles(0),
        levels(NULL), nLevels(0),
        tessellations(NULL), nTessellations(0),
        descendants(NULL), gridID(""),
        gridInputFile("null"), gridOutputFile("null"), gridSoftwareVersion(""),
        gridGenerationDate(""), refCount(0)
    { }
 
    GeoTessGrid(IFStreamAscii& input, const GeoTessOptimizationType* opttype) :
        vertices(NULL), triangles(NULL), levels(NULL), tessellations(NULL),
        descendants(NULL), gridID(""),
        gridInputFile("null"), gridOutputFile("null"), gridSoftwareVersion(""),
        gridGenerationDate(""), refCount(0)
    { loadGridAscii(input); }
 
    GeoTessGrid(IFStreamBinary& input, const GeoTessOptimizationType* opttype) :
        vertices(NULL), triangles(NULL), levels(NULL), tessellations(NULL),
        descendants(NULL), gridID(""),
        gridInputFile("null"), gridOutputFile("null"), gridSoftwareVersion(""),
        gridGenerationDate(""), refCount(0)
    { loadGridBinary(input); }
 
    GeoTessGrid(double** _vertices, int& _nVertices, 
            int** _triangles, int& _nTriangles,
            int** _levels, int& _nLevels,
            int** _tessellations, int& _nTessellations,
            string& _gridID,
            string& _gridInputFile,
            string& _gridOutputFile,
            string& _gridSoftwareVersion,
            string& _gridGenerationDate) :
                vertices(_vertices), nVertices(_nVertices),
                triangles(_triangles), nTriangles(_nTriangles),
                levels(_levels), nLevels(_nLevels),
                tessellations(_tessellations), nTessellations(_nTessellations),
                descendants(NULL),
                gridID(_gridID),
                gridInputFile(_gridInputFile),
                gridOutputFile(_gridOutputFile),
                gridSoftwareVersion(_gridSoftwareVersion),
                gridGenerationDate(_gridGenerationDate),
                refCount(0)
    { initialize(); }
 
    virtual ~GeoTessGrid();
 
    int getReferenceCount() { return refCount; }
 
    void addReference() { ++refCount; }
 
    void removeReference()
    {
        if (isNotReferenced())
        {
            ostringstream os;
            os << endl << "ERROR in GeoTessGrid::removeReference" << endl
                    << "Reference count (" << refCount << ") is already zero." << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 2001);
        }
 
        --refCount;
    }
 
    bool isNotReferenced() { return (refCount == 0); }
 
    void writeGrid(const string& fileName);
 
    void writeGridAscii(const string& fileName);
 
    void writeGridAscii(IFStreamAscii& output);
 
    void writeGridBinary(const string& fileName);
 
    void writeGridBinary(IFStreamBinary& output);
 
    bool isSupportedFormatVersion(int frmtVrsn)
    { return (frmtVrsn == 2); }
 
    // deprecated. GeoTessOptimizationType is no longer used.  Always optimized for speed.
    const GeoTessOptimizationType& getOptimizationType() const { return GeoTessOptimizationType::SPEED; }
 
 
    void testGrid();
 
    int getNLevels() const { return nLevels; }
 
    int const* const* getLevels() const { return levels; }
 
    int const* const* getTessellations() const { return tessellations; }
 
    int getNTriangles() const { return nTriangles; }
 
    vector<Edge*>& getSpokeList(const int& level) const
    { computeSpokeLists(level); return spokeList[level]; }
 
    const vector<vector<Edge*> >& getEdgeList() const { return edgeList; }
 
    const vector<Edge*>& getEdgeList(const int& triangle) const { return edgeList[triangle]; }
 
    void getCenter(const int& triangle, double* center)
    {
        vector<double*> corners;
        corners.push_back(vertices[triangles[triangle][0]]);
        corners.push_back(vertices[triangles[triangle][1]]);
        corners.push_back(vertices[triangles[triangle][2]]);
        GeoTessUtils::center(corners, center);
    }
 
    int delaunay();
 
};
// end class GeoTessGrid
 
}// end namespace geotess
 
#endif  // GEOTESSGRID_OBJECT_H