xypor/index.php

//#     Filename:       SpatialIndex.cpp

//#

//#     The SpatialIndex class is defined here.

//#

//#     Author:         Peter Z. Kunszt based on A. Szalay's code

//#

//#     Date:           October 15, 1998

//#

//#     SPDX-FileCopyrightText: 2000 Peter Z. Kunszt Alex S. Szalay, Aniruddha R. Thakar

//#                     The Johns Hopkins University

//#

//#     Modification History:

//#

//#     Oct 18, 2001 : Dennis C. Dinge -- Replaced ValVec with std::vector

//#     Jul 25, 2002 : Gyorgy Fekete -- Added pointById()

//#


#include "SpatialIndex.h"

#include "SpatialException.h"


#ifdef _WIN32

#include <malloc.h>

#else

#ifdef __APPLE__

#include <sys/malloc.h>

#else

#include <cstdlib>

#endif

#endif


#include <cstdio>

#include <cmath>

// ===========================================================================

//

// Macro definitions for readability

//

// ===========================================================================


#define N(x)      nodes_[(x)]

#define V(x)      vertices_[nodes_[index].v_[(x)]]

#define IV(x)     nodes_[index].v_[(x)]

#define W(x)      vertices_[nodes_[index].w_[(x)]]

#define IW(x)     nodes_[index].w_[(x)]

#define ICHILD(x) nodes_[index].childID_[(x)]


#define IV_(x)     nodes_[index_].v_[(x)]

#define IW_(x)     nodes_[index_].w_[(x)]

#define ICHILD_(x) nodes_[index_].childID_[(x)]

#define IOFFSET    9


// ===========================================================================

//

// Member functions for class SpatialIndex

//

// ===========================================================================


/////////////CONSTRUCTOR//////////////////////////////////

//


SpatialIndex::SpatialIndex(size_t maxlevel, size_t buildlevel)

    : maxlevel_(maxlevel), buildlevel_((buildlevel == 0 || buildlevel > maxlevel) ? maxlevel : buildlevel)

{

    size_t nodes, vertices;


    vMax(&nodes, &vertices);

    layers_.resize(buildlevel_ + 1); // allocate enough space already

    nodes_.resize(nodes + 1);        // allocate space for all nodes

    vertices_.resize(vertices + 1);  // allocate space for all vertices


    N(0).index_ = 0; // initialize invalid node


    // initialize first layer

    layers_[0].level_       = 0;

    layers_[0].nVert_       = 6;

    layers_[0].nNode_       = 8;

    layers_[0].nEdge_       = 12;

    layers_[0].firstIndex_  = 1;

    layers_[0].firstVertex_ = 0;


    // set the first 6 vertices

    float64 v[6][3] = {

        { 0.0L, 0.0L, 1.0L },  // 0

        { 1.0L, 0.0L, 0.0L },  // 1

        { 0.0L, 1.0L, 0.0L },  // 2

        { -1.0L, 0.0L, 0.0L }, // 3

        { 0.0L, -1.0L, 0.0L }, // 4

        { 0.0L, 0.0L, -1.0L }  // 5

    };


    for (int i = 0; i < 6; i++)

        vertices_[i].set(v[i][0], v[i][1], v[i][2]);


    // create the first 8 nodes - index 1 through 8

    index_ = 1;

    newNode(1, 5, 2, 8, 0);  // S0

    newNode(2, 5, 3, 9, 0);  // S1

    newNode(3, 5, 4, 10, 0); // S2

    newNode(4, 5, 1, 11, 0); // S3

    newNode(1, 0, 4, 12, 0); // N0

    newNode(4, 0, 3, 13, 0); // N1

    newNode(3, 0, 2, 14, 0); // N2

    newNode(2, 0, 1, 15, 0); // N3


    //  loop through buildlevel steps, and build the nodes for each layer


    size_t pl    = 0;

    size_t level = buildlevel_;

    while (level-- > 0)

    {

        SpatialEdge edge(*this, pl);

        edge.makeMidPoints();

        makeNewLayer(pl);

        ++pl;

    }


    sortIndex();

}


/////////////NODEVERTEX///////////////////////////////////

// nodeVertex: return the vectors of the vertices, based on the ID

//


void SpatialIndex::nodeVertex(const uint64 id, SpatialVector &v0, SpatialVector &v1, SpatialVector &v2) const

{

    if (buildlevel_ == maxlevel_)

    {

        uint32 idx = (uint32)id - leaves_ + IOFFSET; // -jbb: Fix segfault.  See "idx =" below.

        v0         = vertices_[nodes_[idx].v_[0]];

        v1         = vertices_[nodes_[idx].v_[1]];

        v2         = vertices_[nodes_[idx].v_[2]];

        return;

    }


    // buildlevel < maxlevel

    // get the id of the stored leaf that we are in

    // and get the vertices of the node we want

    uint64 sid = id >> ((maxlevel_ - buildlevel_) * 2);

    uint32 idx = (uint32)(sid - storedleaves_ + IOFFSET);


    v0 = vertices_[nodes_[idx].v_[0]];

    v1 = vertices_[nodes_[idx].v_[1]];

    v2 = vertices_[nodes_[idx].v_[2]];


    // loop through additional levels,

    // pick the correct triangle accordingly, storing the

    // vertices in v1,v2,v3

    for (uint32 i = buildlevel_ + 1; i <= maxlevel_; i++)

    {

        uint64 j         = (id >> ((maxlevel_ - i) * 2)) & 3;

        SpatialVector w0 = v1 + v2;

        w0.normalize();

        SpatialVector w1 = v0 + v2;

        w1.normalize();

        SpatialVector w2 = v1 + v0;

        w2.normalize();


        switch (j)

        {

            case 0:

                v1 = w2;

                v2 = w1;

                break;

            case 1:

                v0 = v1;

                v1 = w0;

                v2 = w2;

                break;

            case 2:

                v0 = v2;

                v1 = w1;

                v2 = w0;

                break;

            case 3:

                v0 = w0;

                v1 = w1;

                v2 = w2;

                break;

        }

    }

}


/////////////MAKENEWLAYER/////////////////////////////////

// makeNewLayer: generate a new layer and the nodes in it

//

void SpatialIndex::makeNewLayer(size_t oldlayer)

{

    uint64 index, id;

    size_t newlayer                = oldlayer + 1;

    layers_[newlayer].level_       = layers_[oldlayer].level_ + 1;

    layers_[newlayer].nVert_       = layers_[oldlayer].nVert_ + layers_[oldlayer].nEdge_;

    layers_[newlayer].nNode_       = 4 * layers_[oldlayer].nNode_;

    layers_[newlayer].nEdge_       = layers_[newlayer].nNode_ + layers_[newlayer].nVert_ - 2;

    layers_[newlayer].firstIndex_  = index_;

    layers_[newlayer].firstVertex_ = layers_[oldlayer].firstVertex_ + layers_[oldlayer].nVert_;


    uint64 ioffset = layers_[oldlayer].firstIndex_;


    for (index = ioffset; index < ioffset + layers_[oldlayer].nNode_; index++)

    {

        id        = N(index).id_ << 2;

        ICHILD(0) = newNode(IV(0), IW(2), IW(1), id++, index);

        ICHILD(1) = newNode(IV(1), IW(0), IW(2), id++, index);

        ICHILD(2) = newNode(IV(2), IW(1), IW(0), id++, index);

        ICHILD(3) = newNode(IW(0), IW(1), IW(2), id, index);

    }

}


/////////////NEWNODE//////////////////////////////////////

// newNode: make a new node

//

uint64 SpatialIndex::newNode(size_t v1, size_t v2, size_t v3, uint64 id, uint64 parent)

{

    IV_(0) = v1; // vertex indices

    IV_(1) = v2;

    IV_(2) = v3;


    IW_(0) = 0; // middle point indices

    IW_(1) = 0;

    IW_(2) = 0;


    ICHILD_(0) = 0; // child indices

    ICHILD_(1) = 0; // index 0 is invalid node.

    ICHILD_(2) = 0;

    ICHILD_(3) = 0;


    N(index_).id_     = id;     // set the id

    N(index_).index_  = index_; // set the index

    N(index_).parent_ = parent; // set the parent


    return index_++;

}


/////////////VMAX/////////////////////////////////////////

// vMax: compute the maximum number of vertices for the

//       polyhedron after buildlevel of subdivisions and

//       the total number of nodes that we store

//       also, calculate the number of leaf nodes that we eventually have.

//

void SpatialIndex::vMax(size_t *nodes, size_t *vertices)

{

    uint64 nv = 6; // initial values

    uint64 ne = 12;

    uint64 nf = 8;

    int32 i   = buildlevel_;

    *nodes    = (size_t)nf;


    while (i-- > 0)

    {

        nv += ne;

        nf *= 4;

        ne = nf + nv - 2;

        *nodes += (size_t)nf;

    }

    *vertices     = (size_t)nv;

    storedleaves_ = nf;


    // calculate number of leaves

    i = maxlevel_ - buildlevel_;

    while (i-- > 0)

        nf *= 4;

    leaves_ = nf;

}


/////////////SORTINDEX////////////////////////////////////

// sortIndex: sort the index so that the first node is the invalid node

//            (index 0), the next 8 nodes are the root nodes

//            and then we put all the leaf nodes in the following block

//            in ascending id-order.

//            All the rest of the nodes is at the end.

void SpatialIndex::sortIndex()

{

    std::vector<QuadNode> oldnodes(nodes_); // create a copy of the node list

    size_t index;

    size_t nonleaf;

    size_t leaf;


#define ON(x) oldnodes[(x)]


    // now refill the nodes_ list according to our sorting.

    for (index = IOFFSET, leaf = IOFFSET, nonleaf = nodes_.size() - 1; index < nodes_.size(); index++)

    {

        if (ON(index).childID_[0] == 0) // childnode

        {

            // set leaf into list

            N(leaf) = ON(index);

            // set parent's pointer to this leaf

            for (size_t i = 0; i < 4; i++)

            {

                if (N(N(leaf).parent_).childID_[i] == index)

                {

                    N(N(leaf).parent_).childID_[i] = leaf;

                    break;

                }

            }

            leaf++;

        }

        else

        {

            // set nonleaf into list from the end

            // set parent of the children already to this

            // index, they come later in the list.

            N(nonleaf)                         = ON(index);

            ON(N(nonleaf).childID_[0]).parent_ = nonleaf;

            ON(N(nonleaf).childID_[1]).parent_ = nonleaf;

            ON(N(nonleaf).childID_[2]).parent_ = nonleaf;

            ON(N(nonleaf).childID_[3]).parent_ = nonleaf;

            // set parent's pointer to this leaf

            for (size_t i = 0; i < 4; i++)

            {

                if (N(N(nonleaf).parent_).childID_[i] == index)

                {

                    N(N(nonleaf).parent_).childID_[i] = nonleaf;

                    break;

                }

            }

            nonleaf--;

        }

    }

}

//////////////////IDBYNAME/////////////////////////////////////////////////

// Translate ascii leaf name to a uint32

//

// The following encoding is used:

//

// The string leaf name has the always the same structure, it begins with

// an N or S, indicating north or south cap and then numbers 0-3 follow

// indicating which child to descend into. So for a depth-5-index we have

// strings like

//                 N012023  S000222  N102302  etc

//

// Each of the numbers correspond to 2 bits of code (00 01 10 11) in the

// uint32. The first two bits are 10 for S and 11 for N. For example

//

//                 N 0 1 2 0 2 3

//                 11000110001011  =  12683 (dec)

//

// The leading bits are always 0.

//

// --- WARNING: This works only up to 15 levels.

//              (we probably never need more than 7)

//


uint64 SpatialIndex::idByName(const char *name)

{

    uint64 out  = 0, i;

    uint32 size = 0;


    if (name == nullptr) // null pointer-name

        throw SpatialFailure("SpatialIndex:idByName:no name given");

    if (name[0] != 'N' && name[0] != 'S') // invalid name

        throw SpatialFailure("SpatialIndex:idByName:invalid name", name);


    size = strlen(name); // determine string length

    // at least size-2 required, don't exceed max

    if (size < 2)

        throw SpatialFailure("SpatialIndex:idByName:invalid name - too short ", name);

    if (size > HTMNAMEMAX)

        throw SpatialFailure("SpatialIndex:idByName:invalid name - too long ", name);


    for (i = size - 1; i > 0; i--) // set bits starting from the end

    {

        if (name[i] > '3' || name[i] < '0') // invalid name

            throw SpatialFailure("SpatialIndex:idByName:invalid name digit ", name);

        out += (uint64(name[i] - '0') << 2 * (size - i - 1));

    }


    i = 2; // set first pair of bits, first bit always set

    if (name[0] == 'N')

        i++; // for north set second bit too

    out += (i << (2 * size - 2));


    /************************

    // This code may be used later for hashing !

    if(size==2)out -= 8;

    else {

      size -= 2;

      uint32 offset = 0, level4 = 8;

      for(i = size; i > 0; i--) { // calculate 4 ^ (level-1), level = size-2

        offset += level4;

        level4 *= 4;

      }

      out -= level4 - offset;

    }

    **************************/

    return out;

}


//////////////////NAMEBYID/////////////////////////////////////////////////

// Translate uint32 to an ascii leaf name

//

// The encoding described above may be decoded again using the following

// procedure:

//

//  * Traverse the uint32 from left to right.

//  * Find the first 'true' bit.

//  * The first pair gives N (11) or S (10).

//  * The subsequent bit-pairs give the numbers 0-3.

//


char *SpatialIndex::nameById(uint64 id, char *name)

{

    uint32 size = 0, i;

#ifdef _WIN32

    uint64 IDHIGHBIT  = 1;

    uint64 IDHIGHBIT2 = 1;

    IDHIGHBIT         = IDHIGHBIT << 63;

    IDHIGHBIT2        = IDHIGHBIT2 << 62;

#endif


    /*************

    // This code might be useful for hashing later !!


    // calculate the level (i.e. 8*4^level) and add it to the id:

    uint32 level=0, level4=8, offset=8;

    while(id >= offset) {

      if(++level > 13) { ok = false; offset = 0; break; }// level too deep

      level4 *= 4;

      offset += level4;

    }

    id += 2 * level4 - offset;

    **************/


    // determine index of first set bit

    for (i = 0; i < IDSIZE; i += 2)

    {

        if ((id << i) & IDHIGHBIT)

            break;

        if ((id << i) & IDHIGHBIT2) // invalid id

            throw SpatialFailure("SpatialIndex:nameById: invalid ID");

    }

    if (id == 0)

        throw SpatialFailure("SpatialIndex:nameById: invalid ID");


    size = (IDSIZE - i) >> 1;

    // allocate characters

    if (!name)

        name = new char[size + 1];


    // fill characters starting with the last one

    for (i = 0; i < size - 1; i++)

        name[size - i - 1] = '0' + char((id >> i * 2) & 3);


    // put in first character

    if ((id >> (size * 2 - 2)) & 1)

    {

        name[0] = 'N';

    }

    else

    {

        name[0] = 'S';

    }

    name[size] = 0; // end string


    return name;

}


//////////////////POINTBYID////////////////////////////////////////////////

// Find a vector for the leaf node given by its ID

//


void SpatialIndex::pointById(SpatialVector &vec, uint64 ID) const

{

    //  uint64 index;

    float64 center_x, center_y, center_z, sum;

    char name[HTMNAMEMAX];


    SpatialVector v0, v1, v2; //


    this->nodeVertex(ID, v0, v1, v2);


    nameById(ID, name);

    /*

        I started to go this way until I discovered nameByID...

        Some docs would be nice for this


        switch(name[1]){

        case '0':

            index = name[0] == 'S' ? 1 : 5;

            break;

        case '1':

            index = name[0] == 'S' ? 2 : 6;

            break;

        case '2':

            index = name[0] == 'S' ? 3 : 7;

            break;

        case '3':

            index = name[0] == 'S' ? 4 : 8;

            break;

        }

    */

    //    cerr << "---------- Point by id: " << name << Qt::endl;

    //  v0.show(); v1.show(); v2.show();

    center_x = v0.x_ + v1.x_ + v2.x_;

    center_y = v0.y_ + v1.y_ + v2.y_;

    center_z = v0.z_ + v1.z_ + v2.z_;

    sum      = center_x * center_x + center_y * center_y + center_z * center_z;

    sum      = sqrt(sum);

    center_x /= sum;

    center_y /= sum;

    center_z /= sum;

    vec.x_ = center_x;

    vec.y_ = center_y;

    vec.z_ = center_z; // I don't want it normalized or radec to be set,

    //  cerr << " - - - - " << Qt::endl;

    //  vec.show();

    //  cerr << "---------- Point by id Retuning" << Qt::endl;

}


//////////////////IDBYPOINT////////////////////////////////////////////////

// Find a leaf node where a vector points to

//


uint64 SpatialIndex::idByPoint(const SpatialVector &v) const

{

    uint64 index;


    // start with the 8 root triangles, find the one which v points to

    for (index = 1; index <= 8; index++)

    {

        if ((V(0) ^ V(1)) * v < -gEpsilon)

            continue;

        if ((V(1) ^ V(2)) * v < -gEpsilon)

            continue;

        if ((V(2) ^ V(0)) * v < -gEpsilon)

            continue;

        break;

    }

    // loop through matching child until leaves are reached

    while (ICHILD(0) != 0)

    {

        uint64 oldindex = index;

        for (size_t i = 0; i < 4; i++)

        {

            index = nodes_[oldindex].childID_[i];

            if ((V(0) ^ V(1)) * v < -gEpsilon)

                continue;

            if ((V(1) ^ V(2)) * v < -gEpsilon)

                continue;

            if ((V(2) ^ V(0)) * v < -gEpsilon)

                continue;

            break;

        }

    }

    // return if we have reached maxlevel

    if (maxlevel_ == buildlevel_)

        return N(index).id_;


    // from now on, continue to build name dynamically.

    // until maxlevel_ levels depth, continue to append the

    // correct index, build the index on the fly.

    char name[HTMNAMEMAX];

    nameById(N(index).id_, name);

    size_t len       = strlen(name);

    SpatialVector v0 = V(0);

    SpatialVector v1 = V(1);

    SpatialVector v2 = V(2);


    size_t level = maxlevel_ - buildlevel_;

    while (level--)

    {

        SpatialVector w0 = v1 + v2;

        w0.normalize();

        SpatialVector w1 = v0 + v2;

        w1.normalize();

        SpatialVector w2 = v1 + v0;

        w2.normalize();


        if (isInside(v, v0, w2, w1))

        {

            name[len++] = '0';

            v1          = w2;

            v2          = w1;

            continue;

        }

        else if (isInside(v, v1, w0, w2))

        {

            name[len++] = '1';

            v0          = v1;

            v1          = w0;

            v2          = w2;

            continue;

        }

        else if (isInside(v, v2, w1, w0))

        {

            name[len++] = '2';

            v0          = v2;

            v1          = w1;

            v2          = w0;

            continue;

        }

        else if (isInside(v, w0, w1, w2))

        {

            name[len++] = '3';

            v0          = w0;

            v1          = w1;

            v2          = w2;

            continue;

        }

    }

    name[len] = '\0';

    return idByName(name);

}


//////////////////ISINSIDE/////////////////////////////////////////////////

// Test whether a vector is inside a triangle. Input triangle has

// to be sorted in a counter-clockwise direction.

//

bool SpatialIndex::isInside(const SpatialVector &v, const SpatialVector &v0, const SpatialVector &v1,

                            const SpatialVector &v2) const

{

    if ((v0 ^ v1) * v < -gEpsilon)

        return false;

    if ((v1 ^ v2) * v < -gEpsilon)

        return false;

    if ((v2 ^ v0) * v < -gEpsilon)

        return false;

    return true;

}

SpatialFailure
SpatialException thrown on operational failure.
Definition SpatialException.h:104

SpatialIndex::idByPoint
uint64 idByPoint(const SpatialVector &vector) const
find a node by giving a vector.
Definition SpatialIndex.cpp:508

SpatialIndex::idByName
static uint64 idByName(const char *)
NodeName conversion to integer ID.
Definition SpatialIndex.cpp:341

SpatialIndex::pointById
void pointById(SpatialVector &vector, uint64 ID) const
find the vector to the centroid of a triangle represented by the ID
Definition SpatialIndex.cpp:457

SpatialIndex::SpatialIndex
SpatialIndex(size_t maxlevel, size_t buildlevel=5)
Constructor.
Definition SpatialIndex.cpp:59

SpatialIndex::nodeVertex
void nodeVertex(const uint64 id, SpatialVector &v1, SpatialVector &v2, SpatialVector &v3) const
return the actual vertex vectors
Definition SpatialIndex.cpp:121

SpatialIndex::nameById
static char * nameById(uint64 ID, char *name=nullptr)
int32 conversion to a string (name of database).
Definition SpatialIndex.cpp:398

SpatialVector
SpatialVector is a 3D vector usually living on the surface of the sphere.
Definition SpatialVector.h:33

SpatialVector::normalize
void normalize()
Normalize vector length to 1.
Definition SpatialVector.cpp:108