← Back to branch summary

~ubuntu-branches/ubuntu/precise/stellarium/precise

« back to all changes in this revision

Viewing changes to src/core/StelGeodesicGrid.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Cédric Delfosse
  • Date: 2009-03-13 20:07:22 UTC
  • mfrom: (1.1.8 upstream) (4.1.2 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090313200722-l66s4zy2s3e8up0s
Tags: 0.10.2-1
New upstream release

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/core/StelGeodesicGrid.cpp
 
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 Vec3d 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 
32
 
 
33
*/
 
34
 
 
35
#include "StelGeodesicGrid.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
StelGeodesicGrid::StelGeodesicGrid(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
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 = NULL;
 
127
}
 
128
 
 
129
void StelGeodesicGrid::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 StelGeodesicGrid::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 StelGeodesicGrid::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 StelGeodesicGrid::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 StelGeodesicGrid::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 StelGeodesicGrid::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 StelGeodesicGrid::searchZone - not found";
 
315
        exit(-1);
 
316
        // shut up the compiler
 
317
        return -1;
 
318
}
 
319
 
 
320
 
 
321
void StelGeodesicGrid::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 StelGeodesicGrid::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* StelGeodesicGrid::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* StelGeodesicGrid::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 StelGeodesicGrid &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[StelGeodesicGrid::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]+StelGeodesicGrid::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
}