1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
/* sphere.c - by David Blythe, SGI */
/* Instead of tessellating a sphere by lines of longitude and latitude
(a technique that over tessellates the poles and under tessellates
the equator of the sphere), tesselate based on regular solids for a
more uniform tesselation.
This approach is arguably better than the gluSphere routine's
approach using slices and stacks (latitude and longitude). -mjk */
#include <stdlib.h>
#include <GL/glut.h>
#include <math.h>
typedef struct {
float x, y, z;
} point;
typedef struct {
point pt[3];
} triangle;
/* six equidistant points lying on the unit sphere */
#define XPLUS { 1, 0, 0 } /* X */
#define XMIN { -1, 0, 0 } /* -X */
#define YPLUS { 0, 1, 0 } /* Y */
#define YMIN { 0, -1, 0 } /* -Y */
#define ZPLUS { 0, 0, 1 } /* Z */
#define ZMIN { 0, 0, -1 } /* -Z */
/* for icosahedron */
#define CZ (0.89442719099991) /* 2/sqrt(5) */
#define SZ (0.44721359549995) /* 1/sqrt(5) */
#define C1 (0.951056516) /* cos(18), */
#define S1 (0.309016994) /* sin(18) */
#define C2 (0.587785252) /* cos(54), */
#define S2 (0.809016994) /* sin(54) */
#define X1 (C1*CZ)
#define Y1 (S1*CZ)
#define X2 (C2*CZ)
#define Y2 (S2*CZ)
#define Ip0 {0., 0., 1.}
#define Ip1 {-X2, -Y2, SZ}
#define Ip2 {X2, -Y2, SZ}
#define Ip3 {X1, Y1, SZ}
#define Ip4 {0, CZ, SZ}
#define Ip5 {-X1, Y1, SZ}
#define Im0 {-X1, -Y1, -SZ}
#define Im1 {0, -CZ, -SZ}
#define Im2 {X1, -Y1, -SZ}
#define Im3 {X2, Y2, -SZ}
#define Im4 {-X2, Y2, -SZ}
#define Im5 {0., 0., -1.}
/* vertices of a unit icosahedron */
static triangle icosahedron[20]= {
/* front pole */
{ {Ip0, Ip1, Ip2}, },
{ {Ip0, Ip5, Ip1}, },
{ {Ip0, Ip4, Ip5}, },
{ {Ip0, Ip3, Ip4}, },
{ {Ip0, Ip2, Ip3}, },
/* mid */
{ {Ip1, Im0, Im1}, },
{ {Im0, Ip1, Ip5}, },
{ {Ip5, Im4, Im0}, },
{ {Im4, Ip5, Ip4}, },
{ {Ip4, Im3, Im4}, },
{ {Im3, Ip4, Ip3}, },
{ {Ip3, Im2, Im3}, },
{ {Im2, Ip3, Ip2}, },
{ {Ip2, Im1, Im2}, },
{ {Im1, Ip2, Ip1}, },
/* back pole */
{ {Im3, Im2, Im5}, },
{ {Im4, Im3, Im5}, },
{ {Im0, Im4, Im5}, },
{ {Im1, Im0, Im5}, },
{ {Im2, Im1, Im5}, },
};
/* normalize point r */
static void
normalize(point *r) {
float mag;
mag = r->x * r->x + r->y * r->y + r->z * r->z;
if (mag != 0.0f) {
mag = 1.0f / sqrt(mag);
r->x *= mag;
r->y *= mag;
r->z *= mag;
}
}
/* linearly interpolate between a & b, by fraction f */
static void
lerp(point *a, point *b, float f, point *r) {
r->x = a->x + f*(b->x-a->x);
r->y = a->y + f*(b->y-a->y);
r->z = a->z + f*(b->z-a->z);
}
void
sphere(int maxlevel) {
int nrows = 1 << maxlevel;
int s;
/* iterate over the 20 sides of the icosahedron */
for(s = 0; s < 20; s++) {
int i;
triangle *t = &icosahedron[s];
for(i = 0; i < nrows; i++) {
/* create a tstrip for each row */
/* number of triangles in this row is number in previous +2 */
/* strip the ith trapezoid block */
point v0, v1, v2, v3, va, vb;
int j;
lerp(&t->pt[1], &t->pt[0], (float)(i+1)/nrows, &v0);
lerp(&t->pt[1], &t->pt[0], (float)i/nrows, &v1);
lerp(&t->pt[1], &t->pt[2], (float)(i+1)/nrows, &v2);
lerp(&t->pt[1], &t->pt[2], (float)i/nrows, &v3);
glBegin(GL_TRIANGLE_STRIP);
#define V(v) { point x; x = v; normalize(&x); glNormal3fv(&x.x); glVertex3fv(&x.x); }
V(v0);
V(v1);
for(j = 0; j < i; j++) {
/* calculate 2 more vertices at a time */
lerp(&v0, &v2, (float)(j+1)/(i+1), &va);
lerp(&v1, &v3, (float)(j+1)/i, &vb);
V(va);
V(vb);
}
V(v2);
#undef V
glEnd();
}
}
}
float *
sphere_tris(int maxlevel) {
int nrows = 1 << maxlevel;
int s, n;
float *buf, *b;
n = 20*(1 << (maxlevel * 2));
b = buf = (float *)malloc(n*3*3*sizeof(float));
/* iterate over the 20 sides of the icosahedron */
for(s = 0; s < 20; s++) {
int i;
triangle *t = &icosahedron[s];
for(i = 0; i < nrows; i++) {
/* create a tstrip for each row */
/* number of triangles in this row is number in previous +2 */
/* strip the ith trapezoid block */
point v0, v1, v2, v3, va, vb, x1, x2;
int j;
lerp(&t->pt[1], &t->pt[0], (float)(i+1)/nrows, &v0);
lerp(&t->pt[1], &t->pt[0], (float)i/nrows, &v1);
lerp(&t->pt[1], &t->pt[2], (float)(i+1)/nrows, &v2);
lerp(&t->pt[1], &t->pt[2], (float)i/nrows, &v3);
#define V(a, c, v) { point x = v; normalize(&a); normalize(&c); normalize(&x); \
b[0] = a.x; b[1] = a.y; b[2] = a.z; \
b[3] = c.x; b[4] = c.y; b[5] = c.z; \
b[6] = x.x; b[7] = x.y; b[8] = x.z; b+=9; }
x1 = v0;
x2 = v1;
for(j = 0; j < i; j++) {
/* calculate 2 more vertices at a time */
lerp(&v0, &v2, (float)(j+1)/(i+1), &va);
lerp(&v1, &v3, (float)(j+1)/i, &vb);
V(x1,x2,va); x1 = x2; x2 = va;
V(vb,x2,x1); x1 = x2; x2 = vb;
}
V(x1, x2, v2);
#undef V
}
}
return buf;
}
|
