← Back to branch summary

~ubuntu-branches/debian/jessie/stellarium/jessie

« back to all changes in this revision

Viewing changes to src/core/GeodesicGrid.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Cédric Delfosse
  • Date: 2009-03-13 20:07:22 UTC
  • mfrom: (1.1.8 upstream)
  • mto: (11.1.1 experimental)
  • mto: This revision was merged to the branch mainline in revision 7.
  • Revision ID: james.westby@ubuntu.com-20090313200722-gbgujsmzsa8a02ty
Import upstream version 0.10.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/core/GeodesicGrid.cpp
1
 
/*
2
 
 
3
 
GeodesicGrid: 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 Vec3d from vecmath.h.
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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
32
 
 
33
 
*/
34
 
 
35
 
#include "GeodesicGrid.hpp"
36
 
 
37
 
#include <QDebug>
38
 
#include <cmath>
39
 
#include <cstdlib>
40
 
 
41
 
static const double icosahedron_G = 0.5*(1.0+sqrt(5.0));
42
 
static const double icosahedron_b = 1.0/sqrt(1.0+icosahedron_G*icosahedron_G);
43
 
static const double icosahedron_a = icosahedron_b*icosahedron_G;
44
 
 
45
 
static const Vec3d icosahedron_corners[12] =
46
 
    {
47
 
        Vec3d( icosahedron_a, -icosahedron_b,            0.0),
48
 
        Vec3d( icosahedron_a,  icosahedron_b,            0.0),
49
 
        Vec3d(-icosahedron_a,  icosahedron_b,            0.0),
50
 
        Vec3d(-icosahedron_a, -icosahedron_b,            0.0),
51
 
        Vec3d(           0.0,  icosahedron_a, -icosahedron_b),
52
 
        Vec3d(           0.0,  icosahedron_a,  icosahedron_b),
53
 
        Vec3d(           0.0, -icosahedron_a,  icosahedron_b),
54
 
        Vec3d(           0.0, -icosahedron_a, -icosahedron_b),
55
 
        Vec3d(-icosahedron_b,            0.0,  icosahedron_a),
56
 
        Vec3d( icosahedron_b,            0.0,  icosahedron_a),
57
 
        Vec3d( icosahedron_b,            0.0, -icosahedron_a),
58
 
        Vec3d(-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
 
GeodesicGrid::GeodesicGrid(const int lev) : maxLevel(lev<0?0:lev), lastMaxSearchlevel(-1)
92
 
{
93
 
        if (maxLevel > 0)
94
 
        {
95
 
                triangles = new Triangle*[maxLevel+1];
96
 
                int 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 = 0;
114
 
        }
115
 
        cacheSearchResult = new GeodesicSearchResult(*this);
116
 
}
117
 
 
118
 
GeodesicGrid::~GeodesicGrid(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 = NULL;
127
 
}
128
 
 
129
 
void GeodesicGrid::getTriangleCorners(int lev,int index,
130
 
                                      Vec3d &h0,
131
 
                                      Vec3d &h1,
132
 
                                      Vec3d &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
 
                switch (index&3)
147
 
                {
148
 
                case 0:
149
 
                                {
150
 
                                    Vec3d c0,c1,c2;
151
 
                                    getTriangleCorners(lev,i,c0,c1,c2);
152
 
                                    h0 = c0;
153
 
                                    h1 = t.e2;
154
 
                                    h2 = t.e1;
155
 
                                }
156
 
                                break;
157
 
                case 1:
158
 
                        {
159
 
                                Vec3d c0,c1,c2;
160
 
                                getTriangleCorners(lev,i,c0,c1,c2);
161
 
                                h0 = t.e2;
162
 
                                h1 = c1;
163
 
                                h2 = t.e0;
164
 
                        }
165
 
                        break;
166
 
                case 2:
167
 
                        {
168
 
                                Vec3d c0,c1,c2;
169
 
                                getTriangleCorners(lev,i,c0,c1,c2);
170
 
                                h0 = t.e1;
171
 
                                h1 = t.e0;
172
 
                                h2 = c2;
173
 
                        }
174
 
                        break;
175
 
                case 3:
176
 
                        h0 = t.e0;
177
 
                        h1 = t.e1;
178
 
                        h2 = t.e2;
179
 
                        break;
180
 
                }
181
 
        }
182
 
}
183
 
 
184
 
int GeodesicGrid::getPartnerTriangle(int lev, int index) const
185
 
{
186
 
        if (lev==0)
187
 
        {
188
 
                Q_ASSERT(index<20);
189
 
                return (index&1) ? index+1 : index-1;
190
 
        }
191
 
        switch(index&7)
192
 
        {
193
 
        case 2:
194
 
        case 6:
195
 
                return index+1;
196
 
        case 3:
197
 
        case 7:
198
 
                return index-1;
199
 
        case 0:
200
 
                return (lev==1) ? index+5 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
201
 
        case 1:
202
 
                return (lev==1) ? index+3 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
203
 
        case 4:
204
 
                return (lev==1) ? index-3 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
205
 
        case 5:
206
 
                return (lev==1) ? index-5 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
207
 
        default:
208
 
                Q_ASSERT(0);
209
 
        }
210
 
        return 0;
211
 
}
212
 
 
213
 
void GeodesicGrid::initTriangle(int lev,int index,
214
 
                                const Vec3d &c0,
215
 
                                const Vec3d &c1,
216
 
                                const Vec3d &c2)
217
 
{
218
 
        Q_ASSERT((c0^c1)*c2 >= 0.0);
219
 
        Triangle &t(triangles[lev][index]);
220
 
        t.e0 = c1+c2;
221
 
        t.e0.normalize();
222
 
        t.e1 = c2+c0;
223
 
        t.e1.normalize();
224
 
        t.e2 = c0+c1;
225
 
        t.e2.normalize();
226
 
        lev++;
227
 
        if (lev < maxLevel)
228
 
        {
229
 
                index *= 4;
230
 
                initTriangle(lev,index+0,c0,t.e2,t.e1);
231
 
                initTriangle(lev,index+1,t.e2,c1,t.e0);
232
 
                initTriangle(lev,index+2,t.e1,t.e0,c2);
233
 
                initTriangle(lev,index+3,t.e0,t.e1,t.e2);
234
 
        }
235
 
}
236
 
 
237
 
 
238
 
void GeodesicGrid::visitTriangles(int maxVisitLevel,
239
 
                                  VisitFunc *func,
240
 
                                  void *context) const
241
 
{
242
 
        if (func && maxVisitLevel >= 0)
243
 
        {
244
 
                if (maxVisitLevel > maxLevel) maxVisitLevel = maxLevel;
245
 
                for (int i=0;i<20;i++)
246
 
                {
247
 
                        const int *const corners = icosahedron_triangles[i].corners;
248
 
                        visitTriangles(0,i,
249
 
                                       icosahedron_corners[corners[0]],
250
 
                                       icosahedron_corners[corners[1]],
251
 
                                       icosahedron_corners[corners[2]],
252
 
                                       maxVisitLevel,func,context);
253
 
                }
254
 
        }
255
 
}
256
 
 
257
 
void GeodesicGrid::visitTriangles(int lev,int index,
258
 
                                  const Vec3d &c0,
259
 
                                  const Vec3d &c1,
260
 
                                  const Vec3d &c2,
261
 
                                  int maxVisitLevel,
262
 
                                  VisitFunc *func,
263
 
                                  void *context) const
264
 
{
265
 
        (*func)(lev,index,c0,c1,c2,context);
266
 
        Triangle &t(triangles[lev][index]);
267
 
        lev++;
268
 
        if (lev <= maxVisitLevel)
269
 
        {
270
 
                index *= 4;
271
 
                visitTriangles(lev,index+0,c0,t.e2,t.e1,maxVisitLevel,func,context);
272
 
                visitTriangles(lev,index+1,t.e2,c1,t.e0,maxVisitLevel,func,context);
273
 
                visitTriangles(lev,index+2,t.e1,t.e0,c2,maxVisitLevel,func,context);
274
 
                visitTriangles(lev,index+3,t.e0,t.e1,t.e2,maxVisitLevel,func,context);
275
 
        }
276
 
}
277
 
 
278
 
 
279
 
int GeodesicGrid::searchZone(const Vec3d &v,int searchLevel) const
280
 
{
281
 
        for (int i=0;i<20;i++)
282
 
        {
283
 
                const int *const corners = icosahedron_triangles[i].corners;
284
 
                const Vec3d &c0(icosahedron_corners[corners[0]]);
285
 
                const Vec3d &c1(icosahedron_corners[corners[1]]);
286
 
                const Vec3d &c2(icosahedron_corners[corners[2]]);
287
 
                if (((c0^c1)*v >= 0.0) && ((c1^c2)*v >= 0.0) && ((c2^c0)*v >= 0.0))
288
 
                {
289
 
                        // v lies inside this icosahedron triangle
290
 
                        for (int lev=0;lev<searchLevel;lev++)
291
 
                        {
292
 
                                Triangle &t(triangles[lev][i]);
293
 
                                i <<= 2;
294
 
                                if ((t.e1^t.e2)*v <= 0.0)
295
 
                                {
296
 
                                        // i += 0;
297
 
                                }
298
 
                                else if ((t.e2^t.e0)*v <= 0.0)
299
 
                                {
300
 
                                        i += 1;
301
 
                                }
302
 
                                else if ((t.e0^t.e1)*v <= 0.0)
303
 
                                {
304
 
                                        i += 2;
305
 
                                }
306
 
                                else
307
 
                                {
308
 
                                        i += 3;
309
 
                                }
310
 
                        }
311
 
                        return i;
312
 
                }
313
 
        }
314
 
        qWarning() << "ERROR GeodesicGrid::searchZone - not found";
315
 
        exit(-1);
316
 
        // shut up the compiler
317
 
        return -1;
318
 
}
319
 
 
320
 
 
321
 
void GeodesicGrid::searchZones(const StelGeom::ConvexS& convex,
322
 
                               int **inside_list,int **border_list,
323
 
                               int maxSearchLevel) const
324
 
{
325
 
        if (maxSearchLevel < 0) maxSearchLevel = 0;
326
 
        else if (maxSearchLevel > maxLevel) maxSearchLevel = maxLevel;
327
 
#if defined __STRICT_ANSI__ || !defined __GNUC__
328
 
        int *halfs_used = new int[convex.size()];
329
 
#else
330
 
        int halfs_used[convex.size()];
331
 
#endif
332
 
        for (int h=0;h<(int)convex.size();h++) {halfs_used[h] = h;}
333
 
#if defined __STRICT_ANSI__ || !defined __GNUC__
334
 
        bool *corner_inside[12];
335
 
        for(int ci=0; ci < 12; ci++) corner_inside[ci]= new bool[convex.size()];
336
 
#else
337
 
        bool corner_inside[12][convex.size()];
338
 
#endif
339
 
        for (size_t h=0;h<convex.size();h++)
340
 
        {
341
 
                const StelGeom::HalfSpace &half_space(convex[h]);
342
 
                for (int i=0;i<12;i++)
343
 
                {
344
 
                        corner_inside[i][h] = half_space.contains(icosahedron_corners[i]);
345
 
                }
346
 
        }
347
 
        for (int i=0;i<20;i++)
348
 
        {
349
 
                searchZones(0,i,
350
 
                            convex,halfs_used,convex.size(),
351
 
                            corner_inside[icosahedron_triangles[i].corners[0]],
352
 
                            corner_inside[icosahedron_triangles[i].corners[1]],
353
 
                            corner_inside[icosahedron_triangles[i].corners[2]],
354
 
                            inside_list,border_list,maxSearchLevel);
355
 
        }
356
 
#if defined __STRICT_ANSI__ || !defined __GNUC__
357
 
        delete[] halfs_used;
358
 
        for(int ci=0; ci < 12; ci++) delete[] corner_inside[ci];
359
 
#endif
360
 
}
361
 
 
362
 
void GeodesicGrid::searchZones(int lev,int index,
363
 
                               const StelGeom::ConvexS& convex,
364
 
                               const int *indexOfUsedHalfSpaces,
365
 
                               const int halfSpacesUsed,
366
 
                               const bool *corner0_inside,
367
 
                               const bool *corner1_inside,
368
 
                               const bool *corner2_inside,
369
 
                               int **inside_list,int **border_list,
370
 
                               const int maxSearchLevel) const
371
 
{
372
 
#if defined __STRICT_ANSI__ || !defined __GNUC__
373
 
        int *halfs_used = new int[halfSpacesUsed];
374
 
#else
375
 
        int halfs_used[halfSpacesUsed];
376
 
#endif
377
 
        int halfs_used_count = 0;
378
 
        for (int h=0;h<halfSpacesUsed;h++)
379
 
        {
380
 
                const int i = indexOfUsedHalfSpaces[h];
381
 
                if (!corner0_inside[i] && !corner1_inside[i] && !corner2_inside[i])
382
 
                {
383
 
                        // totally outside this HalfSpace
384
 
                        return;
385
 
                }
386
 
                else if (corner0_inside[i] && corner1_inside[i] && corner2_inside[i])
387
 
                {
388
 
                        // totally inside this HalfSpace
389
 
                }
390
 
                else
391
 
                {
392
 
                        // on the border of this HalfSpace
393
 
                        halfs_used[halfs_used_count++] = i;
394
 
                }
395
 
        }
396
 
        if (halfs_used_count == 0)
397
 
        {
398
 
                // this triangle(lev,index) lies inside all halfspaces
399
 
                **inside_list = index;
400
 
                (*inside_list)++;
401
 
        }
402
 
        else
403
 
        {
404
 
                (*border_list)--;
405
 
                **border_list = index;
406
 
                if (lev < maxSearchLevel)
407
 
                {
408
 
                        Triangle &t(triangles[lev][index]);
409
 
                        lev++;
410
 
                        index <<= 2;
411
 
                        inside_list++;
412
 
                        border_list++;
413
 
#if defined __STRICT_ANSI__ || !defined __GNUC__
414
 
                        bool *edge0_inside = new bool[convex.size()];
415
 
                        bool *edge1_inside = new bool[convex.size()];
416
 
                        bool *edge2_inside = new bool[convex.size()];
417
 
#else
418
 
                        bool edge0_inside[convex.size()];
419
 
                        bool edge1_inside[convex.size()];
420
 
                        bool edge2_inside[convex.size()];
421
 
#endif
422
 
                        for (int h=0;h<halfs_used_count;h++)
423
 
                        {
424
 
                                const int i = halfs_used[h];
425
 
                                const StelGeom::HalfSpace &half_space(convex[i]);
426
 
                                edge0_inside[i] = half_space.contains(t.e0);
427
 
                                edge1_inside[i] = half_space.contains(t.e1);
428
 
                                edge2_inside[i] = half_space.contains(t.e2);
429
 
                        }
430
 
                        searchZones(lev,index+0,
431
 
                                    convex,halfs_used,halfs_used_count,
432
 
                                    corner0_inside,edge2_inside,edge1_inside,
433
 
                                    inside_list,border_list,maxSearchLevel);
434
 
                        searchZones(lev,index+1,
435
 
                                    convex,halfs_used,halfs_used_count,
436
 
                                    edge2_inside,corner1_inside,edge0_inside,
437
 
                                    inside_list,border_list,maxSearchLevel);
438
 
                        searchZones(lev,index+2,
439
 
                                    convex,halfs_used,halfs_used_count,
440
 
                                    edge1_inside,edge0_inside,corner2_inside,
441
 
                                    inside_list,border_list,maxSearchLevel);
442
 
                        searchZones(lev,index+3,
443
 
                                    convex,halfs_used,halfs_used_count,
444
 
                                    edge0_inside,edge1_inside,edge2_inside,
445
 
                                    inside_list,border_list,maxSearchLevel);
446
 
#if defined __STRICT_ANSI__ || !defined __GNUC__
447
 
                        delete[] edge0_inside;
448
 
                        delete[] edge1_inside;
449
 
                        delete[] edge2_inside;
450
 
#endif
451
 
                }
452
 
        }
453
 
#if defined __STRICT_ANSI__ || !defined __GNUC__
454
 
        delete[] halfs_used;
455
 
#endif
456
 
}
457
 
 
458
 
/*************************************************************************
459
 
 Return a search result matching the given spatial region
460
 
*************************************************************************/
461
 
const GeodesicSearchResult* GeodesicGrid::search(const StelGeom::ConvexS& convex, int maxSearchLevel) const
462
 
{
463
 
        // Try to use the cached version
464
 
        if (maxSearchLevel==lastMaxSearchlevel && convex==lastSearchRegion)
465
 
        {
466
 
                return cacheSearchResult;
467
 
        }
468
 
        // Else recompute it and update cache parameters
469
 
        lastMaxSearchlevel = maxSearchLevel;
470
 
        lastSearchRegion = convex;
471
 
        cacheSearchResult->search(convex, maxSearchLevel);
472
 
        return cacheSearchResult;
473
 
}
474
 
        
475
 
 
476
 
/*************************************************************************
477
 
 Return a search result matching the given spatial region
478
 
*************************************************************************/
479
 
const GeodesicSearchResult* GeodesicGrid::search(const Vec3d &e0,const Vec3d &e1,const Vec3d &e2,const Vec3d &e3,int maxSearchLevel) const
480
 
{
481
 
        StelGeom::ConvexS c(e0, e1, e2, e3);
482
 
        return search(c,maxSearchLevel);
483
 
}
484
 
 
485
 
 
486
 
GeodesicSearchResult::GeodesicSearchResult(const GeodesicGrid &grid)
487
 
                :grid(grid),
488
 
                zones(new int*[grid.getMaxLevel()+1]),
489
 
                inside(new int*[grid.getMaxLevel()+1]),
490
 
                border(new int*[grid.getMaxLevel()+1])
491
 
{
492
 
        for (int i=0;i<=grid.getMaxLevel();i++)
493
 
        {
494
 
                zones[i] = new int[GeodesicGrid::nrOfZones(i)];
495
 
        }
496
 
}
497
 
 
498
 
GeodesicSearchResult::~GeodesicSearchResult(void)
499
 
{
500
 
        for (int i=grid.getMaxLevel();i>=0;i--)
501
 
        {
502
 
                delete[] zones[i];
503
 
        }
504
 
        delete[] border;
505
 
        delete[] inside;
506
 
        delete[] zones;
507
 
}
508
 
 
509
 
void GeodesicSearchResult::search(const StelGeom::ConvexS& convex,
510
 
                                  int maxSearchLevel)
511
 
{
512
 
        for (int i=grid.getMaxLevel();i>=0;i--)
513
 
        {
514
 
                inside[i] = zones[i];
515
 
                border[i] = zones[i]+GeodesicGrid::nrOfZones(i);
516
 
        }
517
 
        grid.searchZones(convex,inside,border,maxSearchLevel);
518
 
}
519
 
 
520
 
void GeodesicSearchInsideIterator::reset(void)
521
 
{
522
 
        level = 0;
523
 
        maxCount = 1<<(maxLevel<<1); // 4^maxLevel
524
 
        indexP = r.zones[0];
525
 
        endP = r.inside[0];
526
 
        index = (*indexP) * maxCount;
527
 
        count = (indexP < endP) ? 0 : maxCount;
528
 
}
529
 
 
530
 
int GeodesicSearchInsideIterator::next(void)
531
 
{
532
 
        if (count < maxCount) return index+(count++);
533
 
        indexP++;
534
 
        if (indexP < endP)
535
 
        {
536
 
                index = (*indexP) * maxCount;
537
 
                count = 1;
538
 
                return index;
539
 
        }
540
 
        while (level < maxLevel)
541
 
        {
542
 
                level++;
543
 
                maxCount >>= 2;
544
 
                indexP = r.zones[level];
545
 
                endP = r.inside[level];
546
 
                if (indexP < endP)
547
 
                {
548
 
                        index = (*indexP) * maxCount;
549
 
                        count = 1;
550
 
                        return index;
551
 
                }
552
 
        }
553
 
        return -1;
554
 
}