15670918640

Committed 16 Jun 2025 02:08AM UTC coverage: 11.775% (-0.2%) from 11.931%

Build # 15670918640

Build Type

push

github

Committed by

alex-w

Commit Message

Updated data

Run Details

14700 of 124846 relevant lines covered (11.77%)

18324.52 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/core/StelGeodesicGrid.cpp

/*
 
StelGeodesicGrid: a library for dividing the sphere into triangle zones
by subdividing the icosahedron
 
Author and Copyright: Johannes Gajdosik, 2006
 
This library requires a simple Vector library,
which may have different copyright and license,
for example Vec3f from VecMath.hpp.
 
In the moment I choose to distribute the library under the GPL,
later I may choose to additionally distribute it under a more
relaxed license like the LGPL. If you want to have the library
under another license, please ask me.
 
 
 
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
 
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
 
You should have received a copy of the GNU General Public License
along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA  02110-1335, USA.
 
*/

#include "StelGeodesicGrid.hpp"

#include <QDebug>
#include <cmath>
#include <cstdlib>

static const float icosahedron_G = 0.5f*(1.0f+std::sqrt(5.0f));
static const float icosahedron_b = 1.0f/std::sqrt(1.0f+icosahedron_G*icosahedron_G);
static const float icosahedron_a = icosahedron_b*icosahedron_G;

static const Vec3f icosahedron_corners[12] =
    {
                Vec3f( icosahedron_a, -icosahedron_b,            0.0),
                Vec3f( icosahedron_a,  icosahedron_b,            0.0),
                Vec3f(-icosahedron_a,  icosahedron_b,            0.0),
                Vec3f(-icosahedron_a, -icosahedron_b,            0.0),
                Vec3f(           0.0,  icosahedron_a, -icosahedron_b),
                Vec3f(           0.0,  icosahedron_a,  icosahedron_b),
                Vec3f(           0.0, -icosahedron_a,  icosahedron_b),
                Vec3f(           0.0, -icosahedron_a, -icosahedron_b),
                Vec3f(-icosahedron_b,            0.0,  icosahedron_a),
                Vec3f( icosahedron_b,            0.0,  icosahedron_a),
                Vec3f( icosahedron_b,            0.0, -icosahedron_a),
                Vec3f(-icosahedron_b,            0.0, -icosahedron_a)
    };

struct TopLevelTriangle
{
        int corners[3];   // index der Ecken
};


static const TopLevelTriangle icosahedron_triangles[20] =
    {
        {{ 1, 0,10}}, //  1
        {{ 0, 1, 9}}, //  0
        {{ 0, 9, 6}}, // 12
        {{ 9, 8, 6}}, //  9
        {{ 0, 7,10}}, // 16
        {{ 6, 7, 0}}, //  6
        {{ 7, 6, 3}}, //  7
        {{ 6, 8, 3}}, // 14
        {{11,10, 7}}, // 11
        {{ 7, 3,11}}, // 18
        {{ 3, 2,11}}, //  3
        {{ 2, 3, 8}}, //  2
        {{10,11, 4}}, // 10
        {{ 2, 4,11}}, // 19
        {{ 5, 4, 2}}, //  5
        {{ 2, 8, 5}}, // 15
        {{ 4, 1,10}}, // 17
        {{ 4, 5, 1}}, //  4
        {{ 5, 9, 1}}, // 13
        {{ 8, 9, 5}}  //  8
    };

StelGeodesicGrid::StelGeodesicGrid(const int lev) : maxLevel(lev<0?0:lev), lastMaxSearchlevel(-1)
{
        if (maxLevel > 0)
        {
                triangles = new Triangle*[static_cast<size_t>(maxLevel)+1];
                uint nr_of_triangles = 20;
                for (int i=0;i<maxLevel;i++)
                {
                        triangles[i] = new Triangle[nr_of_triangles];
                        nr_of_triangles *= 4;
                }
                for (int i=0;i<20;i++)
                {
                        const int *const corners = icosahedron_triangles[i].corners;
                        initTriangle(0,i,
                                     icosahedron_corners[corners[0]],
                                     icosahedron_corners[corners[1]],
                                     icosahedron_corners[corners[2]]);
                }
        }
        else
        {
                triangles = Q_NULLPTR;
        }
        cacheSearchResult = new GeodesicSearchResult(*this);
}

StelGeodesicGrid::~StelGeodesicGrid(void)
{
        if (maxLevel > 0)
        {
                for (int i=maxLevel-1;i>=0;i--) delete[] triangles[i];
                delete[] triangles;
        }
        delete cacheSearchResult;
        cacheSearchResult = Q_NULLPTR;
}

void StelGeodesicGrid::getTriangleCorners(int lev,int index,
                                                                          Vec3f &h0,
                                                                          Vec3f &h1,
                                                                          Vec3f &h2) const
{
        if (lev <= 0)
        {
                const int *const corners = icosahedron_triangles[index].corners;
                h0 = icosahedron_corners[corners[0]];
                h1 = icosahedron_corners[corners[1]];
                h2 = icosahedron_corners[corners[2]];
        }
        else
        {
                lev--;
                const int i = index>>2;
                Triangle &t(triangles[lev][i]);
                Vec3f c0,c1,c2;
                switch (index&3)
                {
                        case 0:
                                getTriangleCorners(lev,i,c0,c1,c2);
                                h0 = c0;
                                h1 = t.e2;
                                h2 = t.e1;
                                break;
                        case 1:
                                getTriangleCorners(lev,i,c0,c1,c2);
                                h0 = t.e2;
                                h1 = c1;
                                h2 = t.e0;
                                break;
                        case 2:
                                getTriangleCorners(lev,i,c0,c1,c2);
                                h0 = t.e1;
                                h1 = t.e0;
                                h2 = c2;
                                break;
                        case 3:
                                h0 = t.e0;
                                h1 = t.e1;
                                h2 = t.e2;
                                break;
                }
        }
}

int StelGeodesicGrid::getPartnerTriangle(int lev, int index) const
{
        if (lev==0)
        {
                Q_ASSERT(index<20);
                return (index&1) ? index+1 : index-1;
        }
        switch(index&7)
        {
                case 2:
                case 6:
                        return index+1;
                case 3:
                case 7:
                        return index-1;
                case 0:
                        return (lev==1) ? index+5 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
                case 1:
                        return (lev==1) ? index+3 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
                case 4:
                        return (lev==1) ? index-3 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
                case 5:
                        return (lev==1) ? index-5 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
                default:
                        Q_ASSERT(0);
        }
        return 0;
}

void StelGeodesicGrid::initTriangle(int lev,int index,
                                                                const Vec3f &c0,
                                                                const Vec3f &c1,
                                                                const Vec3f &c2)
{
        Q_ASSERT((c0^c1)*c2 >= 0.0f);
        Triangle &t(triangles[lev][index]);
        t.e0 = c1+c2;
        t.e0.normalize();
        t.e1 = c2+c0;
        t.e1.normalize();
        t.e2 = c0+c1;
        t.e2.normalize();
        lev++;
        if (lev < maxLevel)
        {
                index *= 4;
                initTriangle(lev,index+0,c0,t.e2,t.e1);
                initTriangle(lev,index+1,t.e2,c1,t.e0);
                initTriangle(lev,index+2,t.e1,t.e0,c2);
                initTriangle(lev,index+3,t.e0,t.e1,t.e2);
        }
}


void StelGeodesicGrid::visitTriangles(int maxVisitLevel,
                                  VisitFunc *func,
                                  void *context) const
{
        if (func && maxVisitLevel >= 0)
        {
                if (maxVisitLevel > maxLevel) maxVisitLevel = maxLevel;
                for (int i=0;i<20;i++)
                {
                        const int *const corners = icosahedron_triangles[i].corners;
                        visitTriangles(0,i,
                                       icosahedron_corners[corners[0]],
                                       icosahedron_corners[corners[1]],
                                       icosahedron_corners[corners[2]],
                                       maxVisitLevel,func,context);
                }
        }
}

void StelGeodesicGrid::visitTriangles(int lev,int index,
                                                                  const Vec3f &c0,
                                                                  const Vec3f &c1,
                                                                  const Vec3f &c2,
                                  int maxVisitLevel,
                                  VisitFunc *func,
                                  void *context) const
{
        (*func)(lev,index,c0,c1,c2,context);
        Triangle &t(triangles[lev][index]);
        lev++;
        if (lev <= maxVisitLevel)
        {
                index *= 4;
                visitTriangles(lev,index+0,c0,t.e2,t.e1,maxVisitLevel,func,context);
                visitTriangles(lev,index+1,t.e2,c1,t.e0,maxVisitLevel,func,context);
                visitTriangles(lev,index+2,t.e1,t.e0,c2,maxVisitLevel,func,context);
                visitTriangles(lev,index+3,t.e0,t.e1,t.e2,maxVisitLevel,func,context);
        }
}


int StelGeodesicGrid::getZoneNumberForPoint(const Vec3f &v,int searchLevel) const
{
        for (int i=0;i<20;i++)
        {
                const int *const corners = icosahedron_triangles[i].corners;
                const Vec3f &c0(icosahedron_corners[corners[0]]);
                const Vec3f &c1(icosahedron_corners[corners[1]]);
                const Vec3f &c2(icosahedron_corners[corners[2]]);
                if (((c0^c1)*v >= 0.0f) && ((c1^c2)*v >= 0.0f) && ((c2^c0)*v >= 0.0f))
                {
                        // v lies inside this icosahedron triangle
                        for (int lev=0;lev<searchLevel;lev++)
                        {
                                Triangle &t(triangles[lev][i]);
                                i <<= 2;
                                if ((t.e1^t.e2)*v <= 0.0f)
                                {
                                        // i += 0;
                                }
                                else if ((t.e2^t.e0)*v <= 0.0f)
                                {
                                        i += 1;
                                }
                                else if ((t.e0^t.e1)*v <= 0.0f)
                                {
                                        i += 2;
                                }
                                else
                                {
                                        i += 3;
                                }
                        }
                        return i;
                }
        }
        qWarning() << "ERROR StelGeodesicGrid::searchZone - not found";
        exit(-1);
        // shut up the compiler
        return -1;
}


// First iteration on the icosahedron base triangles
void StelGeodesicGrid::searchZones(const QVector<SphericalCap>& convex,
                               int **inside_list,int **border_list,
                               int maxSearchLevel) const
{
        if (maxSearchLevel < 0) maxSearchLevel = 0;
        else if (maxSearchLevel > maxLevel) maxSearchLevel = maxLevel;
#if defined __STRICT_ANSI__ || !defined __GNUC__
        int *halfs_used = new int[static_cast<size_t>(convex.size())];
#else
        int halfs_used[static_cast<size_t>(convex.size())];
#endif
        for (int h=0;h<static_cast<int>(convex.size());h++) {halfs_used[h] = h;}
#if defined __STRICT_ANSI__ || !defined __GNUC__
        bool *corner_inside[12];
        for(int ci=0; ci < 12; ci++) corner_inside[ci]= new bool[static_cast<size_t>(convex.size())];
#else
        bool corner_inside[12][static_cast<size_t>(convex.size())];
#endif
        for (int h=0;h<convex.size();h++)
        {
                const SphericalCap& half_space(convex.at(h));
                for (int i=0;i<12;i++)
                {
                        corner_inside[i][h] = half_space.contains(icosahedron_corners[i]);
                }
        }
        for (int i=0;i<20;i++)
        {
                searchZones(0,i,
                            convex,halfs_used,convex.size(),
                            corner_inside[icosahedron_triangles[i].corners[0]],
                            corner_inside[icosahedron_triangles[i].corners[1]],
                            corner_inside[icosahedron_triangles[i].corners[2]],
                            inside_list,border_list,maxSearchLevel);
        }
#if defined __STRICT_ANSI__ || !defined __GNUC__
        delete[] halfs_used;
        for(int ci=0; ci < 12; ci++) delete[] corner_inside[ci];
#endif
}

void StelGeodesicGrid::searchZones(int lev,int index,
                                                                   const QVector<SphericalCap>&convex,
                               const int *indexOfUsedSphericalCaps,
                               const int halfSpacesUsed,
                               const bool *corner0_inside,
                               const bool *corner1_inside,
                               const bool *corner2_inside,
                               int **inside_list,int **border_list,
                               const int maxSearchLevel) const
{
#if defined __STRICT_ANSI__ || !defined __GNUC__
        int *halfs_used = new int[static_cast<size_t>(halfSpacesUsed)];
#else
        int halfs_used[static_cast<size_t>(halfSpacesUsed)];
#endif
        int halfs_used_count = 0;
        for (int h=0;h<halfSpacesUsed;h++)
        {
                const int i = indexOfUsedSphericalCaps[h];
                if (!corner0_inside[i] && !corner1_inside[i] && !corner2_inside[i])
                {
                        // totally outside this SphericalCap
                        goto end;
                }
                else if (corner0_inside[i] && corner1_inside[i] && corner2_inside[i])
                {
                        // totally inside this SphericalCap
                }
                else
                {
                        // on the border of this SphericalCap
                        halfs_used[halfs_used_count++] = i;
                }
        }
        if (halfs_used_count == 0)
        {
                // this triangle(lev,index) lies inside all halfspaces
                **inside_list = index;
                (*inside_list)++;
        }
        else
        {
                (*border_list)--;
                **border_list = index;
                if (lev < maxSearchLevel)
                {
                        const Triangle &t(triangles[lev][index]);
                        lev++;
                        index <<= 2;
                        inside_list++;
                        border_list++;
#if defined __STRICT_ANSI__ || !defined __GNUC__
                        bool *edge0_inside = new bool[static_cast<size_t>(convex.size())];
                        bool *edge1_inside = new bool[static_cast<size_t>(convex.size())];
                        bool *edge2_inside = new bool[static_cast<size_t>(convex.size())];
#else
                        bool edge0_inside[static_cast<size_t>(convex.size())];
                        bool edge1_inside[static_cast<size_t>(convex.size())];
                        bool edge2_inside[static_cast<size_t>(convex.size())];
#endif
                        for (int h=0;h<halfs_used_count;h++)
                        {
                                const int i = halfs_used[h];
                                const SphericalCap& half_space(convex.at(i));
                                edge0_inside[i] = half_space.contains(t.e0);
                                edge1_inside[i] = half_space.contains(t.e1);
                                edge2_inside[i] = half_space.contains(t.e2);
                        }
                        searchZones(lev,index+0,
                                    convex,halfs_used,halfs_used_count,
                                    corner0_inside,edge2_inside,edge1_inside,
                                    inside_list,border_list,maxSearchLevel);
                        searchZones(lev,index+1,
                                    convex,halfs_used,halfs_used_count,
                                    edge2_inside,corner1_inside,edge0_inside,
                                    inside_list,border_list,maxSearchLevel);
                        searchZones(lev,index+2,
                                    convex,halfs_used,halfs_used_count,
                                    edge1_inside,edge0_inside,corner2_inside,
                                    inside_list,border_list,maxSearchLevel);
                        searchZones(lev,index+3,
                                    convex,halfs_used,halfs_used_count,
                                    edge0_inside,edge1_inside,edge2_inside,
                                    inside_list,border_list,maxSearchLevel);
#if defined __STRICT_ANSI__ || !defined __GNUC__
                        delete[] edge0_inside;
                        delete[] edge1_inside;
                        delete[] edge2_inside;
#endif
                }
        }
end:
#if defined __STRICT_ANSI__ || !defined __GNUC__
        delete[] halfs_used;
#endif
        return;
}

/*************************************************************************
 Return a search result matching the given spatial region
*************************************************************************/
const GeodesicSearchResult* StelGeodesicGrid::search(const QVector<SphericalCap>& convex, int maxSearchLevel) const
{
        // Try to use the cached version
        if (maxSearchLevel==lastMaxSearchlevel && convex==lastSearchRegion)
        {
                return cacheSearchResult;
        }
        // Else recompute it and update cache parameters
        lastMaxSearchlevel = maxSearchLevel;
        lastSearchRegion = convex;
        cacheSearchResult->search(convex, maxSearchLevel);
        return cacheSearchResult;
}


GeodesicSearchResult::GeodesicSearchResult(const StelGeodesicGrid &grid)
                :grid(grid),
                zones( new int*[static_cast<size_t>(grid.getMaxLevel())+1]),
                inside(new int*[static_cast<size_t>(grid.getMaxLevel())+1]),
                border(new int*[static_cast<size_t>(grid.getMaxLevel())+1])
{
        for (int i=0;i<=grid.getMaxLevel();i++)
        {
                zones[i] = new int[static_cast<size_t>(StelGeodesicGrid::nrOfZones(i))];
        }
}

GeodesicSearchResult::~GeodesicSearchResult(void)
{
        for (int i=grid.getMaxLevel();i>=0;i--)
        {
                delete[] zones[i];
        }
        delete[] border;
        delete[] inside;
        delete[] zones;
}

void GeodesicSearchResult::print() const
{
        qDebug() << "GeodesicSearchResult: Grid max level:" << grid.getMaxLevel() << "Grid nr of zones:" << grid.getNrOfZones() <<
                    "zones have " << grid.getMaxLevel()+1 << "entries";
        for (int i=0;i<=grid.getMaxLevel();i++)
        {
                qDebug() << "     Zone " << i << ": " << static_cast<size_t>(StelGeodesicGrid::nrOfZones(i)) << "entries";
        }
}

void GeodesicSearchResult::search(const QVector<SphericalCap>& convex, int maxSearchLevel)
{
        for (int i=grid.getMaxLevel();i>=0;i--)
        {
                inside[i] = zones[i];
                border[i] = zones[i]+StelGeodesicGrid::nrOfZones(i);
        }
        grid.searchZones(convex,inside,border,maxSearchLevel);
}

void GeodesicSearchInsideIterator::reset(void)
{
        level = 0;
        maxCount = 1<<(maxLevel<<1); // 4^maxLevel
        indexP = r.zones[0];
        endP = r.inside[0];
        index = (*indexP) * maxCount;
        count = (indexP < endP) ? 0 : maxCount;
}

int GeodesicSearchInsideIterator::next(void)
{
        if (count < maxCount) return index+(count++);
        indexP++;
        if (indexP < endP)
        {
                index = (*indexP) * maxCount;
                count = 1;
                return index;
        }
        while (level < maxLevel)
        {
                level++;
                maxCount >>= 2;
                indexP = r.zones[level];
                endP = r.inside[level];
                if (indexP < endP)
                {
                        index = (*indexP) * maxCount;
                        count = 1;
                        return index;
                }
        }
        return -1;
}

1	/*
2
3	StelGeodesicGrid: a library for dividing the sphere into triangle zones
4	by subdividing the icosahedron
5
6	Author and Copyright: Johannes Gajdosik, 2006
7
8	This library requires a simple Vector library,
9	which may have different copyright and license,
10	for example Vec3f from VecMath.hpp.
11
12	In the moment I choose to distribute the library under the GPL,
13	later I may choose to additionally distribute it under a more
14	relaxed license like the LGPL. If you want to have the library
15	under another license, please ask me.
16
17
18
19	This library is free software; you can redistribute it and/or
20	modify it under the terms of the GNU General Public License
21	as published by the Free Software Foundation; either version 2
22	of the License, or (at your option) any later version.
23
24	This library is distributed in the hope that it will be useful,
25	but WITHOUT ANY WARRANTY; without even the implied warranty of
26	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27	GNU General Public License for more details.
28
29	You should have received a copy of the GNU General Public License
30	along with this library; if not, write to the Free Software
31	Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
32
33	*/
34
35	#include "StelGeodesicGrid.hpp"
36
37	#include <QDebug>
38	#include <cmath>
39	#include <cstdlib>
40
41	static const float icosahedron_G = 0.5f*(1.0f+std::sqrt(5.0f));
42	static const float icosahedron_b = 1.0f/std::sqrt(1.0f+icosahedron_G*icosahedron_G);
43	static const float icosahedron_a = icosahedron_b*icosahedron_G;
44
45	static const Vec3f icosahedron_corners[12] =
46	{
47	Vec3f( icosahedron_a, -icosahedron_b, 0.0),
48	Vec3f( icosahedron_a, icosahedron_b, 0.0),
49	Vec3f(-icosahedron_a, icosahedron_b, 0.0),
50	Vec3f(-icosahedron_a, -icosahedron_b, 0.0),
51	Vec3f( 0.0, icosahedron_a, -icosahedron_b),
52	Vec3f( 0.0, icosahedron_a, icosahedron_b),
53	Vec3f( 0.0, -icosahedron_a, icosahedron_b),
54	Vec3f( 0.0, -icosahedron_a, -icosahedron_b),
55	Vec3f(-icosahedron_b, 0.0, icosahedron_a),
56	Vec3f( icosahedron_b, 0.0, icosahedron_a),
57	Vec3f( icosahedron_b, 0.0, -icosahedron_a),
58	Vec3f(-icosahedron_b, 0.0, -icosahedron_a)
59	};
60
61	struct TopLevelTriangle
62	{
63	int corners[3]; // index der Ecken
64	};
65
66
67	static const TopLevelTriangle icosahedron_triangles[20] =
68	{
69	{{ 1, 0,10}}, // 1
70	{{ 0, 1, 9}}, // 0
71	{{ 0, 9, 6}}, // 12
72	{{ 9, 8, 6}}, // 9
73	{{ 0, 7,10}}, // 16
74	{{ 6, 7, 0}}, // 6
75	{{ 7, 6, 3}}, // 7
76	{{ 6, 8, 3}}, // 14
77	{{11,10, 7}}, // 11
78	{{ 7, 3,11}}, // 18
79	{{ 3, 2,11}}, // 3
80	{{ 2, 3, 8}}, // 2
81	{{10,11, 4}}, // 10
82	{{ 2, 4,11}}, // 19
83	{{ 5, 4, 2}}, // 5
84	{{ 2, 8, 5}}, // 15
85	{{ 4, 1,10}}, // 17
86	{{ 4, 5, 1}}, // 4
87	{{ 5, 9, 1}}, // 13
88	{{ 8, 9, 5}} // 8
89	};
90
91	StelGeodesicGrid::StelGeodesicGrid(const int lev) : maxLevel(lev<0?0:lev), lastMaxSearchlevel(-1)	×
92	{
93	if (maxLevel > 0)	×
94	{
95	triangles = new Triangle*[static_cast<size_t>(maxLevel)+1];	×
96	uint nr_of_triangles = 20;	×
97	for (int i=0;i<maxLevel;i++)	×
98	{
99	triangles[i] = new Triangle[nr_of_triangles];	×
100	nr_of_triangles *= 4;	×
101	}
102	for (int i=0;i<20;i++)	×
103	{
104	const int *const corners = icosahedron_triangles[i].corners;	×
105	initTriangle(0,i,	×
106	icosahedron_corners[corners[0]],	×
107	icosahedron_corners[corners[1]],	×
108	icosahedron_corners[corners[2]]);	×
109	}
110	}
111	else
112	{
113	triangles = Q_NULLPTR;	×
114	}
115	cacheSearchResult = new GeodesicSearchResult(*this);	×
116	}	×
117
118	StelGeodesicGrid::~StelGeodesicGrid(void)	×
119	{
120	if (maxLevel > 0)	×
121	{
122	for (int i=maxLevel-1;i>=0;i--) delete[] triangles[i];	×
123	delete[] triangles;	×
124	}
125	delete cacheSearchResult;	×
126	cacheSearchResult = Q_NULLPTR;	×
127	}	×
128
129	void StelGeodesicGrid::getTriangleCorners(int lev,int index,	×
130	Vec3f &h0,
131	Vec3f &h1,
132	Vec3f &h2) const
133	{
134	if (lev <= 0)	×
135	{
136	const int *const corners = icosahedron_triangles[index].corners;	×
137	h0 = icosahedron_corners[corners[0]];	×
138	h1 = icosahedron_corners[corners[1]];	×
139	h2 = icosahedron_corners[corners[2]];	×
140	}
141	else
142	{
143	lev--;	×
144	const int i = index>>2;	×
145	Triangle &t(triangles[lev][i]);	×
146	Vec3f c0,c1,c2;	×
147	switch (index&3)	×
148	{
149	case 0:	×
150	getTriangleCorners(lev,i,c0,c1,c2);	×
151	h0 = c0;	×
152	h1 = t.e2;	×
153	h2 = t.e1;	×
154	break;	×
155	case 1:	×
156	getTriangleCorners(lev,i,c0,c1,c2);	×
157	h0 = t.e2;	×
158	h1 = c1;	×
159	h2 = t.e0;	×
160	break;	×
161	case 2:	×
162	getTriangleCorners(lev,i,c0,c1,c2);	×
163	h0 = t.e1;	×
164	h1 = t.e0;	×
165	h2 = c2;	×
166	break;	×
167	case 3:	×
168	h0 = t.e0;	×
169	h1 = t.e1;	×
170	h2 = t.e2;	×
171	break;	×
172	}
173	}
174	}	×
175
176	int StelGeodesicGrid::getPartnerTriangle(int lev, int index) const	×
177	{
178	if (lev==0)	×
179	{
180	Q_ASSERT(index<20);	×
181	return (index&1) ? index+1 : index-1;	×
182	}
183	switch(index&7)	×
184	{
185	case 2:	×
186	case 6:
187	return index+1;	×
188	case 3:	×
189	case 7:
190	return index-1;	×
191	case 0:	×
192	return (lev==1) ? index+5 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;	×
193	case 1:	×
194	return (lev==1) ? index+3 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;	×
195	case 4:	×
196	return (lev==1) ? index-3 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;	×
197	case 5:	×
198	return (lev==1) ? index-5 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;	×
199	default:	×
200	Q_ASSERT(0);	×
201	}
202	return 0;
203	}
204
205	void StelGeodesicGrid::initTriangle(int lev,int index,	×
206	const Vec3f &c0,
207	const Vec3f &c1,
208	const Vec3f &c2)
209	{
210	Q_ASSERT((c0^c1)*c2 >= 0.0f);	×
211	Triangle &t(triangles[lev][index]);	×
212	t.e0 = c1+c2;	×
213	t.e0.normalize();	×
214	t.e1 = c2+c0;	×
215	t.e1.normalize();	×
216	t.e2 = c0+c1;	×
217	t.e2.normalize();	×
218	lev++;	×
219	if (lev < maxLevel)	×
220	{
221	index *= 4;	×
222	initTriangle(lev,index+0,c0,t.e2,t.e1);	×
223	initTriangle(lev,index+1,t.e2,c1,t.e0);	×
224	initTriangle(lev,index+2,t.e1,t.e0,c2);	×
225	initTriangle(lev,index+3,t.e0,t.e1,t.e2);	×
226	}
227	}	×
228
229
230	void StelGeodesicGrid::visitTriangles(int maxVisitLevel,	×
231	VisitFunc *func,
232	void *context) const
233	{
234	if (func && maxVisitLevel >= 0)	×
235	{
236	if (maxVisitLevel > maxLevel) maxVisitLevel = maxLevel;	×
237	for (int i=0;i<20;i++)	×
238	{
239	const int *const corners = icosahedron_triangles[i].corners;	×
240	visitTriangles(0,i,	×
241	icosahedron_corners[corners[0]],	×
242	icosahedron_corners[corners[1]],	×
243	icosahedron_corners[corners[2]],	×
244	maxVisitLevel,func,context);
245	}
246	}
247	}	×
248
249	void StelGeodesicGrid::visitTriangles(int lev,int index,	×
250	const Vec3f &c0,
251	const Vec3f &c1,
252	const Vec3f &c2,
253	int maxVisitLevel,
254	VisitFunc *func,
255	void *context) const
256	{
257	(*func)(lev,index,c0,c1,c2,context);	×
258	Triangle &t(triangles[lev][index]);	×
259	lev++;	×
260	if (lev <= maxVisitLevel)	×
261	{
262	index *= 4;	×
263	visitTriangles(lev,index+0,c0,t.e2,t.e1,maxVisitLevel,func,context);	×
264	visitTriangles(lev,index+1,t.e2,c1,t.e0,maxVisitLevel,func,context);	×
265	visitTriangles(lev,index+2,t.e1,t.e0,c2,maxVisitLevel,func,context);	×
266	visitTriangles(lev,index+3,t.e0,t.e1,t.e2,maxVisitLevel,func,context);	×
267	}
268	}	×
269
270
271	int StelGeodesicGrid::getZoneNumberForPoint(const Vec3f &v,int searchLevel) const	×
272	{
273	for (int i=0;i<20;i++)	×
274	{
275	const int *const corners = icosahedron_triangles[i].corners;	×
276	const Vec3f &c0(icosahedron_corners[corners[0]]);	×
277	const Vec3f &c1(icosahedron_corners[corners[1]]);	×
278	const Vec3f &c2(icosahedron_corners[corners[2]]);	×
279	if (((c0^c1)v >= 0.0f) && ((c1^c2)v >= 0.0f) && ((c2^c0)*v >= 0.0f))	×
280	{
281	// v lies inside this icosahedron triangle
282	for (int lev=0;lev<searchLevel;lev++)	×
283	{
284	Triangle &t(triangles[lev][i]);	×
285	i <<= 2;	×
286	if ((t.e1^t.e2)*v <= 0.0f)	×
287	{
288	// i += 0;
289	}
290	else if ((t.e2^t.e0)*v <= 0.0f)	×
291	{
292	i += 1;	×
293	}
294	else if ((t.e0^t.e1)*v <= 0.0f)	×
295	{
296	i += 2;	×
297	}
298	else
299	{
300	i += 3;	×
301	}
302	}
303	return i;	×
304	}
305	}
306	qWarning() << "ERROR StelGeodesicGrid::searchZone - not found";	×
307	exit(-1);	×
308	// shut up the compiler
309	return -1;
310	}
311
312
313	// First iteration on the icosahedron base triangles
314	void StelGeodesicGrid::searchZones(const QVector<SphericalCap>& convex,	×
315	int inside_list,int border_list,
316	int maxSearchLevel) const
317	{
318	if (maxSearchLevel < 0) maxSearchLevel = 0;	×
319	else if (maxSearchLevel > maxLevel) maxSearchLevel = maxLevel;	×
320	#if defined __STRICT_ANSI__ \|\| !defined __GNUC__
321	int *halfs_used = new int[static_cast<size_t>(convex.size())];	×
322	#else
323	int halfs_used[static_cast<size_t>(convex.size())];
324	#endif
325	for (int h=0;h<static_cast<int>(convex.size());h++) {halfs_used[h] = h;}	×
326	#if defined __STRICT_ANSI__ \|\| !defined __GNUC__
327	bool *corner_inside[12];
328	for(int ci=0; ci < 12; ci++) corner_inside[ci]= new bool[static_cast<size_t>(convex.size())];	×
329	#else
330	bool corner_inside[12][static_cast<size_t>(convex.size())];
331	#endif
332	for (int h=0;h<convex.size();h++)	×
333	{
334	const SphericalCap& half_space(convex.at(h));	×
335	for (int i=0;i<12;i++)	×
336	{
337	corner_inside[i][h] = half_space.contains(icosahedron_corners[i]);	×
338	}
339	}
340	for (int i=0;i<20;i++)	×
341	{
342	searchZones(0,i,	×
343	convex,halfs_used,convex.size(),	×
344	corner_inside[icosahedron_triangles[i].corners[0]],	×
345	corner_inside[icosahedron_triangles[i].corners[1]],	×
346	corner_inside[icosahedron_triangles[i].corners[2]],	×
347	inside_list,border_list,maxSearchLevel);
348	}
349	#if defined __STRICT_ANSI__ \|\| !defined __GNUC__
350	delete[] halfs_used;	×
351	for(int ci=0; ci < 12; ci++) delete[] corner_inside[ci];	×
352	#endif
353	}	×
354
355	void StelGeodesicGrid::searchZones(int lev,int index,	×
356	const QVector<SphericalCap>&convex,
357	const int *indexOfUsedSphericalCaps,
358	const int halfSpacesUsed,
359	const bool *corner0_inside,
360	const bool *corner1_inside,
361	const bool *corner2_inside,
362	int inside_list,int border_list,
363	const int maxSearchLevel) const
364	{
365	#if defined __STRICT_ANSI__ \|\| !defined __GNUC__
366	int *halfs_used = new int[static_cast<size_t>(halfSpacesUsed)];	×
367	#else
368	int halfs_used[static_cast<size_t>(halfSpacesUsed)];
369	#endif
370	int halfs_used_count = 0;	×
371	for (int h=0;h<halfSpacesUsed;h++)	×
372	{
373	const int i = indexOfUsedSphericalCaps[h];	×
374	if (!corner0_inside[i] && !corner1_inside[i] && !corner2_inside[i])	×
375	{
376	// totally outside this SphericalCap
377	goto end;	×
378	}
379	else if (corner0_inside[i] && corner1_inside[i] && corner2_inside[i])	×
380	{
381	// totally inside this SphericalCap
382	}
383	else
384	{
385	// on the border of this SphericalCap
386	halfs_used[halfs_used_count++] = i;	×
387	}
388	}
389	if (halfs_used_count == 0)	×
390	{
391	// this triangle(lev,index) lies inside all halfspaces
392	**inside_list = index;	×
393	(*inside_list)++;	×
394	}
395	else
396	{
397	(*border_list)--;	×
398	**border_list = index;	×
399	if (lev < maxSearchLevel)	×
400	{
401	const Triangle &t(triangles[lev][index]);	×
402	lev++;	×
403	index <<= 2;	×
404	inside_list++;	×
405	border_list++;	×
406	#if defined __STRICT_ANSI__ \|\| !defined __GNUC__
407	bool *edge0_inside = new bool[static_cast<size_t>(convex.size())];	×
408	bool *edge1_inside = new bool[static_cast<size_t>(convex.size())];	×
409	bool *edge2_inside = new bool[static_cast<size_t>(convex.size())];	×
410	#else
411	bool edge0_inside[static_cast<size_t>(convex.size())];
412	bool edge1_inside[static_cast<size_t>(convex.size())];
413	bool edge2_inside[static_cast<size_t>(convex.size())];
414	#endif
415	for (int h=0;h<halfs_used_count;h++)	×
416	{
417	const int i = halfs_used[h];	×
418	const SphericalCap& half_space(convex.at(i));	×
419	edge0_inside[i] = half_space.contains(t.e0);	×
420	edge1_inside[i] = half_space.contains(t.e1);	×
421	edge2_inside[i] = half_space.contains(t.e2);	×
422	}
423	searchZones(lev,index+0,	×
424	convex,halfs_used,halfs_used_count,
425	corner0_inside,edge2_inside,edge1_inside,
426	inside_list,border_list,maxSearchLevel);
427	searchZones(lev,index+1,	×
428	convex,halfs_used,halfs_used_count,
429	edge2_inside,corner1_inside,edge0_inside,
430	inside_list,border_list,maxSearchLevel);
431	searchZones(lev,index+2,	×
432	convex,halfs_used,halfs_used_count,
433	edge1_inside,edge0_inside,corner2_inside,
434	inside_list,border_list,maxSearchLevel);
435	searchZones(lev,index+3,	×
436	convex,halfs_used,halfs_used_count,
437	edge0_inside,edge1_inside,edge2_inside,
438	inside_list,border_list,maxSearchLevel);
439	#if defined __STRICT_ANSI__ \|\| !defined __GNUC__
440	delete[] edge0_inside;	×
441	delete[] edge1_inside;	×
442	delete[] edge2_inside;	×
443	#endif
444	}
445	}
446	end:	×
447	#if defined __STRICT_ANSI__ \|\| !defined __GNUC__
448	delete[] halfs_used;	×
449	#endif
450	return;	×
451	}
452
453	/*************************************************************************
454	Return a search result matching the given spatial region
455	*************************************************************************/
456	const GeodesicSearchResult* StelGeodesicGrid::search(const QVector<SphericalCap>& convex, int maxSearchLevel) const	×
457	{
458	// Try to use the cached version
459	if (maxSearchLevel==lastMaxSearchlevel && convex==lastSearchRegion)	×
460	{
461	return cacheSearchResult;	×
462	}
463	// Else recompute it and update cache parameters
464	lastMaxSearchlevel = maxSearchLevel;	×
465	lastSearchRegion = convex;	×
466	cacheSearchResult->search(convex, maxSearchLevel);	×
467	return cacheSearchResult;	×
468	}
469
470
471	GeodesicSearchResult::GeodesicSearchResult(const StelGeodesicGrid &grid)	×
472	:grid(grid),	×
473	zones( new int*[static_cast<size_t>(grid.getMaxLevel())+1]),	×
474	inside(new int*[static_cast<size_t>(grid.getMaxLevel())+1]),	×
475	border(new int*[static_cast<size_t>(grid.getMaxLevel())+1])	×
476	{
477	for (int i=0;i<=grid.getMaxLevel();i++)	×
478	{
479	zones[i] = new int[static_cast<size_t>(StelGeodesicGrid::nrOfZones(i))];	×
480	}
481	}	×
482
483	GeodesicSearchResult::~GeodesicSearchResult(void)	×
484	{
485	for (int i=grid.getMaxLevel();i>=0;i--)	×
486	{
487	delete[] zones[i];	×
488	}
489	delete[] border;	×
490	delete[] inside;	×
491	delete[] zones;	×
492	}	×
493
494	void GeodesicSearchResult::print() const	×
495	{
496	qDebug() << "GeodesicSearchResult: Grid max level:" << grid.getMaxLevel() << "Grid nr of zones:" << grid.getNrOfZones() <<	×
497	"zones have " << grid.getMaxLevel()+1 << "entries";	×
498	for (int i=0;i<=grid.getMaxLevel();i++)	×
499	{
500	qDebug() << " Zone " << i << ": " << static_cast<size_t>(StelGeodesicGrid::nrOfZones(i)) << "entries";	×
501	}
502	}	×
503
504	void GeodesicSearchResult::search(const QVector<SphericalCap>& convex, int maxSearchLevel)	×
505	{
506	for (int i=grid.getMaxLevel();i>=0;i--)	×
507	{
508	inside[i] = zones[i];	×
509	border[i] = zones[i]+StelGeodesicGrid::nrOfZones(i);	×
510	}
511	grid.searchZones(convex,inside,border,maxSearchLevel);	×
512	}	×
513
514	void GeodesicSearchInsideIterator::reset(void)	×
515	{
516	level = 0;	×
517	maxCount = 1<<(maxLevel<<1); // 4^maxLevel	×
518	indexP = r.zones[0];	×
519	endP = r.inside[0];	×
520	index = (indexP) maxCount;	×
521	count = (indexP < endP) ? 0 : maxCount;	×
522	}	×
523
524	int GeodesicSearchInsideIterator::next(void)	×
525	{
526	if (count < maxCount) return index+(count++);	×
527	indexP++;	×
528	if (indexP < endP)	×
529	{
530	index = (indexP) maxCount;	×
531	count = 1;	×
532	return index;	×
533	}
534	while (level < maxLevel)	×
535	{
536	level++;	×
537	maxCount >>= 2;	×
538	indexP = r.zones[level];	×
539	endP = r.inside[level];	×
540	if (indexP < endP)	×
541	{
542	index = (indexP) maxCount;	×
543	count = 1;	×
544	return index;	×
545	}
546	}
547	return -1;	×
548	}

Stellarium / stellarium / 15670918640

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous