1 /* Copyright (c) 2015  Gerald Knizia
2  *
3  * This file is part of the IboView program (see: http://www.iboview.org)
4  *
5  * IboView is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, version 3.
8  *
9  * IboView is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
16  *
17  * Please see IboView documentation in README.txt for:
18  * -- A list of included external software and their licenses. The included
19  *    external software's copyright is not touched by this agreement.
20  * -- Notes on re-distribution and contributions to/further development of
21  *    the IboView software
22  */
23 
24 #include "Iv.h"
25 #include <algorithm>
26 #include <GL/glew.h>
27 #include "IvMesh.h"
28 #include "CxPodArray.h"
29 #include "CxVec3.h"
30 
31 
32 template<class FVertex>
TIndexedTriangleList(FVec3d * Pos,FVec3d * Norm,size_t nVertices,unsigned * Indices,size_t nIndices,FVertex const & RefVertex)33 TIndexedTriangleList<FVertex>::TIndexedTriangleList(FVec3d *Pos, FVec3d *Norm, size_t nVertices, unsigned *Indices, size_t nIndices, FVertex const &RefVertex)
34 {
35    // translate to target location and convert to GL format
36    Vertices.reserve(nVertices);
37    for (unsigned iVert = 0; iVert < nVertices; ++ iVert) {
38       FVertex
39          vOut(RefVertex);
40       FVec3d const
41          &vIn = Pos[iVert],
42          &vNorm = Norm[iVert];
43       vOut.vNorm = vec3f(vNorm[0], vNorm[1], vNorm[2]);
44       vOut.vPos  = vec3f(vIn[0], vIn[1], vIn[2]);
45       Vertices.push_back(vOut);
46    }
47 
48    assert(nIndices % 3 == 0);
49    Triangles.reserve(nIndices);
50    for (unsigned iTri = 0; iTri < nIndices/3; ++ iTri) {
51       Triangles.push_back(vec3ui(Indices[3*iTri+0], Indices[3*iTri+1], Indices[3*iTri+2]));
52    }
53 }
54 
55 
56 template<class FVertex>
Append(TIndexedTriangleList const & Other)57 void TIndexedTriangleList<FVertex>::Append(TIndexedTriangleList const &Other)
58 {
59    unsigned
60       iVertexOffset = Vertices.size(),
61       iTriangle0 = Triangles.size();
62    Vertices.insert(Vertices.end(), Other.Vertices.begin(), Other.Vertices.end());
63    Triangles.insert(Triangles.end(), Other.Triangles.begin(), Other.Triangles.end());
64    // fix up triangle indices of new triangles.
65    for (unsigned i = iTriangle0; i < Triangles.size(); ++ i)
66       for (uint iVx = 0; iVx < 3; ++ iVx)
67          Triangles[i][iVx] += iVertexOffset;
68 }
69 
70 
71 namespace SubdivSphere {
72    using ct::TArray;
73    typedef TArray<FTriangle1>
74       FTriangleArray;
75 
76    // see Icosahedron.py
77    static const double fIcosahedronCoordinates[12][3] = {
78       {0.85065080835203999, 0.0, 0.52573111211913359},
79       {0.85065080835203999, 0.0, -0.52573111211913359},
80       {-0.85065080835203999, 0.0, 0.52573111211913359},
81       {-0.85065080835203999, 0.0, -0.52573111211913359},
82       {0.0, 0.52573111211913359, 0.85065080835203999},
83       {0.0, -0.52573111211913359, 0.85065080835203999},
84       {0.0, 0.52573111211913359, -0.85065080835203999},
85       {0.0, -0.52573111211913359, -0.85065080835203999},
86       {0.52573111211913359, 0.85065080835203999, 0.0},
87       {-0.52573111211913359, 0.85065080835203999, 0.0},
88       {0.52573111211913359, -0.85065080835203999, 0.0},
89       {-0.52573111211913359, -0.85065080835203999, 0.0}
90    };
91    static const unsigned iIcosahedronTriangles[20][3] = {{0,1,8},{0,4,5},{0,5,10},{0,8,4},{0,10,1},
92       {1,6,8},{1,7,6},{1,10,7},{2,3,11},{2,4,9},{2,5,4},{2,9,3},{2,11,5},{3,6,7},{3,7,11},{3,9,6},
93       {4,8,9},{5,11,10},{6,9,8},{7,10,11}};
94 
95 
96    // make a non-indexed triangle list out of a given list of positions and indices.
MakeInitialMesh(FTriangleArray & Out,double const (* vCoords)[3],unsigned const (* iTriangles)[3],unsigned nTriangles)97    void MakeInitialMesh(FTriangleArray &Out, double const (*vCoords)[3], unsigned const (*iTriangles)[3], unsigned nTriangles) {
98       Out.reserve(nTriangles);
99       for (unsigned iTri = 0; iTri < nTriangles; ++ iTri) {
100          unsigned const
101             *ii = &iTriangles[iTri][0];
102          FTriangle1 t;
103          t.v[0] = FVec3d(vCoords[ii[0]]);
104          t.v[1] = FVec3d(vCoords[ii[1]]);
105          t.v[2] = FVec3d(vCoords[ii[2]]);
106          Out.push_back(t);
107       }
108    }
109 
110    // subdivide the given triangles by bisecting their edges and projecting
111    // the new points onto the unit sphere.
112    //
113    // (note: there are some valid subdivisions which cannot be accessed with this:
114    //  technically we could divide the edges into any number of equal length segments,
115    //  and project those back.)
BinarilySubdivideUnitSphere(FTriangleArray & Out,FTriangleArray const & In)116    void BinarilySubdivideUnitSphere(FTriangleArray &Out, FTriangleArray const &In)
117    {
118       Out.clear();
119       Out.reserve(In.size() * 4);
120       for (unsigned iTri = 0; iTri < In.size(); ++ iTri) {
121          //            0
122          //
123          //
124          //      5          3
125          //
126          //
127          //  2        4        1
128          FVec3d v[6];
129          for (unsigned i = 0; i < 3; ++i)
130             v[i] = In[iTri].v[i];
131          for (unsigned i = 3; i < 6; ++i) {
132             v[i] = .5 * (v[i-3] + v[(i-2)%3]); // e.g., new vertex 3 goes between old vertices 0 and 1.
133             // project new coordinate to unit sphere.
134             v[i] /= ct::Length(v[i]);
135          }
136          Out.push_back(FTriangle1(v[5], v[0], v[3]));
137          Out.push_back(FTriangle1(v[1], v[4], v[3]));
138          Out.push_back(FTriangle1(v[3], v[4], v[5]));
139          Out.push_back(FTriangle1(v[5], v[4], v[2]));
140       }
141    }
142 
143    // make non-indexed triangle list representing a unit sphere by recursively
144    // sub-dividing an icosahedron.
MakeSubdivUnitSphere(FTriangleArray & Out,unsigned nSubdivLevel)145    void MakeSubdivUnitSphere(FTriangleArray &Out, unsigned nSubdivLevel)
146    {
147       // take an icosahedron and divide it along the edges, each time normalizing
148       // the triangles back onto the sphere.
149       FTriangleArray
150          T0, T1;
151       MakeInitialMesh(T0, fIcosahedronCoordinates, iIcosahedronTriangles, 20);
152       for (unsigned iSubdivLevel = 0; iSubdivLevel < nSubdivLevel; ++ iSubdivLevel) {
153          // subdivide T0 into T1.
154          BinarilySubdivideUnitSphere(T1, T0);
155          // exchange them.
156          T1.swap(T0);
157       }
158       Out.swap(T0);
159    }
160 
161    // define some auxiliary structures for joining equal vertices on different triangles.
162    // (to turn a direct triangle list into an indexed triangle list)
163    static const FVec3d RandomPlane(1.2344591,0.0823435367,-.323451231);
164    static const double fEpsilon = 1e-8;
165 
166    struct FVertexInfo {
167       FVec3d vPos;
168       unsigned iOrig;
169       double fSortPred;
FVertexInfoSubdivSphere::FVertexInfo170       FVertexInfo(FVec3d vPos_, unsigned iOrig_)
171          : vPos(vPos_), iOrig(iOrig_), fSortPred(ct::Dot(RandomPlane, vPos_))
172       {}
173 
174       // induce a weak ordering to simplify searching for equal vertices.
operator <SubdivSphere::FVertexInfo175       bool operator < (FVertexInfo const &other) const {
176          // sort by distance to a random plane.
177          return this->fSortPred < other.fSortPred;
178       }
179 
operator ==SubdivSphere::FVertexInfo180       bool operator == (FVertexInfo const &other) const {
181          // are both are sufficiently close to be considered equal?
182          return DistSq(this->vPos, other.vPos) < fEpsilon*fEpsilon;
183          // note: this test is not compatible with operator <.
184       }
185    };
186 
187 
188    // find unique positions in PosIn and return them in PosOut.
189    // Indices is an array of length PosIn.size() such that PosOut[Indices[i]] == PosIn[i].
FindUniquePositions(TArray<FVec3d> & PosOut,TArray<unsigned> & Indices,TArray<FVec3d> const & PosIn)190    void FindUniquePositions(TArray<FVec3d> &PosOut, TArray<unsigned> &Indices, TArray<FVec3d> const &PosIn)
191    {
192       TArray<FVertexInfo>
193          vs;
194       vs.reserve(PosIn.size());
195       for (unsigned i = 0; i < PosIn.size(); ++ i)
196          vs.push_back(FVertexInfo(PosIn[i], i));
197       std::sort(vs.begin(), vs.end());
198 
199       PosOut.clear();
200       PosOut.reserve(PosIn.size());
201       Indices.resize(PosIn.size());
202       unsigned
203          iNext;
204       for (unsigned i = 0; i < vs.size(); i = iNext) {
205          iNext = i + 1;
206          while (iNext != vs.size() && vs[i] == vs[iNext])
207             iNext += 1;
208          // ^- note: that is not quite right. Vertices which are too far apart
209          // according to our sorting criterion cannot be equal, but vertices
210          // which are equal need non necessarily end up rigt next to each other
211          // in the metric. we should search through *all* vertices which lie close
212          // to the current one according to the sorting metric.
213          for (unsigned k = i; k != iNext; ++ k)
214             Indices[vs[k].iOrig] = PosOut.size();
215          PosOut.push_back(vs[i].vPos);
216       }
217    }
218 
219 
TranslateAndScale(ct::TArray<FVec3d> & Pos,double fScale,FVec3d vDisplacement)220    void TranslateAndScale(ct::TArray<FVec3d> &Pos, double fScale, FVec3d vDisplacement) {
221       for (unsigned i = 0; i < Pos.size(); ++ i)
222          Pos[i] = fScale * Pos[i] + vDisplacement;
223    }
224 
225 } // namespace SubdivSphere
226 
227 
MakeIndexedTriangles(ct::TArray<FVec3d> & PosOut,ct::TArray<unsigned> & Indices,ct::TArray<FTriangle1> const & TrianglesIn)228 void MakeIndexedTriangles(ct::TArray<FVec3d> &PosOut, ct::TArray<unsigned> &Indices, ct::TArray<FTriangle1> const &TrianglesIn)
229 {
230    using namespace SubdivSphere;
231 
232    // flatten the triangle array to turn it into a list of position vectors.
233    TArray<FVec3d>
234       Triangles;
235    Triangles.reserve(3*TrianglesIn.size());
236    for (unsigned iTri = 0; iTri < TrianglesIn.size(); ++ iTri)
237       for (unsigned iv = 0; iv != 3; ++ iv)
238          Triangles.push_back(TrianglesIn[iTri].v[iv]);
239    // make indexed triangles.
240    TArray<unsigned>
241       VertexIndices;
242    FindUniquePositions(PosOut, VertexIndices, Triangles);
243    Indices.clear();
244    Indices.reserve(3*TrianglesIn.size());
245    for (unsigned iTri = 0; iTri < TrianglesIn.size(); ++ iTri)
246       for (unsigned iv = 0; iv != 3; ++ iv)
247          Indices.push_back(VertexIndices[3*iTri + iv]);
248 }
249 
250 
MakeSubdivSphere(vec3f vPos,float fRadius,uint nSubdivLevel,FBaseVertex const & RefVertex)251 FIndexedTriangleListPtr MakeSubdivSphere(vec3f vPos, float fRadius, uint nSubdivLevel, FBaseVertex const &RefVertex)
252 {
253    using namespace SubdivSphere;
254 
255    // make vertex positions and triangle indices of unit sphere
256    ct::TArray<FVec3d> Pos;
257    ct::TArray<unsigned> Indices;
258    {
259       FTriangleArray
260          T0;
261       MakeSubdivUnitSphere(T0, nSubdivLevel);
262       MakeIndexedTriangles(Pos, Indices, T0);
263    }
264 
265    // translate to target location and convert to GL format
266    TranslateAndScale(Pos, (double)fRadius, FVec3d(vPos[0], vPos[1], vPos[2]));
267    return FIndexedTriangleListPtr(new FIndexedTriangleList(&Pos[0], &Pos[0], Pos.size(), &Indices[0], Indices.size(), RefVertex));
268 }
269 
MakeCylinder(float fBaseRadius,float fTopRadius,float fHeight,uint nSlices,FBaseVertex const & RefVertex)270 FIndexedTriangleListPtr MakeCylinder(float fBaseRadius, float fTopRadius, float fHeight, uint nSlices, FBaseVertex const &RefVertex)
271 {
272    using namespace SubdivSphere;
273 
274    // make vertex positions and triangle indices of cylinder
275    ct::TArray<FVec3d> Pos, Norm;
276    ct::TArray<unsigned> Indices;
277 
278    unsigned n = nSlices;
279    // base and top circle
280    Pos.resize(2*n);
281    Norm.resize(2*n);
282    // a circle with nSlices has nSlices edges, and each of those leads to
283    // one face composed of two triangles.
284    Indices.resize(3*2*n);
285 
286    //       \         t
287    //     ny \     /<-->|
288    //         \   /     |
289    //       <->\ /      |
290    //        nx /       |
291    //          /        |
292    //         /<------->|
293    //              b
294 
295    double
296       nx = fHeight,
297       ny = -(fTopRadius - fBaseRadius),
298       nf = 1./std::sqrt(nx*nx + ny*ny);
299    nx *= nf;
300    ny *= nf;
301 
302    for (unsigned iAngle = 0; iAngle < n; ++ iAngle) {
303       double
304          fAngle = 2 * M_PI * iAngle/float(n),
305          cx = std::cos(fAngle),
306          cy = std::sin(fAngle);
307       Pos[iAngle] = FVec3d(fBaseRadius*cx, fBaseRadius*cy, 0.);
308       Pos[iAngle+n] = FVec3d(fTopRadius*cx, fTopRadius*cy, fHeight);
309       Norm[iAngle] = FVec3d(nx*cx, nx*cy, ny);
310       Norm[iAngle+n] = FVec3d(nx*cx, nx*cy, ny);
311 //       Norm[iAngle+n] = FVec3d(cx, cy, 0.); // <- to give a bent look?
312 
313       {
314          //    iAng-----iAng-1
315          //    |    /
316          //    |   /
317          //    iAng+n   iAng+n-1
318          Indices[6*iAngle+0] = iAngle;
319          Indices[6*iAngle+1] = (iAngle + 1)%n;
320          Indices[6*iAngle+2] = iAngle + n;
321          Indices[6*iAngle+3] = iAngle + n;
322          Indices[6*iAngle+4] = (iAngle + 1)%n;
323          Indices[6*iAngle+5] = (iAngle + 1)%n + n;
324       }
325    }
326    return FIndexedTriangleListPtr(new FIndexedTriangleList(&Pos[0], &Norm[0], Pos.size(), &Indices[0], Indices.size(), RefVertex));
327 }
328 
MakeIcosahedron(float fRadius,FBaseVertex const & RefVertex)329 FIndexedTriangleListPtr MakeIcosahedron(float fRadius, FBaseVertex const &RefVertex)
330 {
331    using namespace SubdivSphere;
332    // we don't bother with indexing it up here. we just create three vertices
333    // and three indices per face.
334    ct::TArray<FVec3d> Pos, Norm;
335    ct::TArray<unsigned> Indices;
336 
337    unsigned nFaces = 20;
338    Pos.resize(3*nFaces);
339    Norm.resize(3*nFaces);
340    Indices.reserve(3*nFaces);
341 
342    for (size_t iFace = 0; iFace < 20; ++ iFace) {
343       // we don't bother with indexing it up. we just create three vertices and three indices
344       // per face.
345       FVec3d
346          &v0 = Pos[3*iFace],
347          &v1 = Pos[3*iFace+1],
348          &v2 = Pos[3*iFace+2];
349 
350       unsigned const
351          *pFaceIndices = iIcosahedronTriangles[iFace];
352       for (size_t ixyz = 0; ixyz < 3; ++ ixyz) {
353          v0[ixyz] = fRadius * fIcosahedronCoordinates[pFaceIndices[0]][ixyz];
354          v1[ixyz] = fRadius * fIcosahedronCoordinates[pFaceIndices[1]][ixyz];
355          v2[ixyz] = fRadius * fIcosahedronCoordinates[pFaceIndices[2]][ixyz];
356       }
357       FVec3d
358          vFaceNormal = Normalized(Cross(v1-v0,v2-v0));
359       for (size_t iVert = 0; iVert < 3; ++ iVert) {
360          Norm[3*iFace + iVert] = vFaceNormal;
361          Indices.push_back(3*iFace + iVert);
362       }
363    }
364    assert(Pos.size() == 60 && Norm.size() == 60 && Indices.size() == 60);
365    return FIndexedTriangleListPtr(new FIndexedTriangleList(&Pos[0], &Norm[0], Pos.size(), &Indices[0], Indices.size(), RefVertex));
366 }
367 
368 
369 
370 template struct TIndexedTriangleList<FBaseVertex>;
371