1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23 /****************************************************************************
24 History
25
26 $Log: gen_normal.h,v $
27 ****************************************************************************/
28
29 #ifndef __VCG_GEN_NORMAL
30 #define __VCG_GEN_NORMAL
31
32 #include <algorithm>
33
34 namespace vcg {
35
36 template <class ScalarType>
37 class GenNormal
38 {
39 public:
40 typedef Point3<ScalarType> Point3x;
41
Random(int vn,std::vector<Point3<ScalarType>> & NN)42 static void Random(int vn, std::vector<Point3<ScalarType > > &NN)
43 {
44 NN.clear();
45 while(NN.size()<vn)
46 {
47 Point3x pp(((float)rand())/RAND_MAX,
48 ((float)rand())/RAND_MAX,
49 ((float)rand())/RAND_MAX);
50 pp=pp*2.0-Point3x(1,1,1);
51 if(pp.SquaredNorm()<1)
52 {
53 Normalize(pp);
54 NN.push_back(pp);
55 }
56 }
57 }
58
59
FibonacciPt(int i,int n)60 static Point3x FibonacciPt(int i, int n)
61 {
62 const ScalarType Phi = ScalarType(std::sqrt(ScalarType(5))*0.5 + 0.5);
63 const ScalarType phi = 2.0*M_PI* (i/Phi - floor(i/Phi));
64 ScalarType cosTheta = 1.0 - (2*i + 1.0)/ScalarType(n);
65 ScalarType sinTheta = 1 - cosTheta*cosTheta;
66 sinTheta = std::sqrt(std::min(ScalarType(1),std::max(ScalarType(0),sinTheta)));
67 return Point3x(
68 cos(phi)*sinTheta,
69 sin(phi)*sinTheta,
70 cosTheta);
71 }
72
73 // Implementation of the Spherical Fibonacci Point Sets
74 // according to the description of
75 // Spherical Fibonacci Mapping
76 // Benjamin Keinert, Matthias Innmann, Michael Sanger, Marc Stamminger
77 // TOG 2015
Fibonacci(int n,std::vector<Point3x> & NN)78 static void Fibonacci(int n, std::vector<Point3x > &NN)
79 {
80 NN.resize(n);
81 for(int i=0;i<n;++i)
82 NN[i]=FibonacciPt(i,n);
83 }
84
85 static void UniformCone(int vn, std::vector<Point3<ScalarType > > &NN, ScalarType AngleRad, Point3x dir=Point3x(0,1,0))
86 {
87 std::vector<Point3<ScalarType > > NNT;
88 NN.clear();
89 // per prima cosa si calcola il volume della spherical cap di angolo AngleRad
90 ScalarType Height= 1.0 - cos(AngleRad); // height is measured from top...
91 // Surface is the one of the tangent cylinder
92 ScalarType CapArea = 2.0*M_PI*Height;
93 ScalarType Ratio = CapArea / (4.0*M_PI );
94
95 printf("----------AngleRad %f Angledeg %f ratio %f vn %i vn2 %i \n",AngleRad,math::ToDeg(AngleRad),Ratio,vn,int(vn/Ratio));
96 Fibonacci(vn/Ratio,NNT);
97 printf("asked %i got %i (expecting %i instead of %i)\n", int(vn/Ratio), NNT.size(), int(NNT.size()*Ratio), vn);
98 typename std::vector<Point3<ScalarType> >::iterator vi;
99
100 ScalarType cosAngle = cos(AngleRad);
101 for(vi=NNT.begin();vi!=NNT.end();++vi)
102 {
103 if(dir.dot(*vi) >= cosAngle) NN.push_back(*vi);
104 }
105 }
106
107 // This is an Implementation of the Dave Rusin’s Disco Ball algorithm
108 // You can spread the points as follows:
109 // Put N+1 points on the meridian from north to south poles, equally spaced.
110 // If you swing this meridian around the sphere, you'll sweep out the entire
111 // surface; in the process, each of the points will sweep out a circle. You
112 // can show that the ith point will sweep out a circle of radius sin(pi i/N).
113 // If you space points equally far apart on this circle, keeping the
114 // displacement roughly the same as on that original meridian, you'll be
115 // able to fit about 2N sin(pi i/N) points here. This process will put points
116 // pretty evenly spaced on the sphere; the number of such points is about
117 // 2+ 2N*Sum(i=1 to N-1) sin(pi i/N).
118 // The closed form of this summation
119 // 2.0 - ( (2.0*N * sin (M_PI/N))/(cos(M_PI/N) - 1.0));
DiscoBall(int vn,std::vector<Point3<ScalarType>> & NN)120 static void DiscoBall(int vn, std::vector<Point3<ScalarType > > &NN)
121 {
122 // Guess the right N
123 ScalarType N=0;
124
125 for(N=1;N<vn;++N)
126 {
127 ScalarType expectedPoints = 2.0 - ( (2.0*N * sin (M_PI/N))/(cos(M_PI/N) - 1.0));
128 if(expectedPoints >= vn) break;
129 }
130
131 ScalarType VerticalAngle = M_PI / N;
132 NN.push_back(Point3<ScalarType>(0,0,1.0));
133 for (int i =1; i<N; ++i)
134 {
135 // Z is the north/south axis
136 ScalarType HorizRadius = sin(i*VerticalAngle);
137 ScalarType CircleLength = 2.0 * M_PI * HorizRadius;
138 ScalarType Z = cos(i*VerticalAngle);
139 ScalarType PointNumPerCircle = floor( CircleLength / VerticalAngle);
140 ScalarType HorizontalAngle = 2.0*M_PI/PointNumPerCircle;
141 for(ScalarType j=0;j<PointNumPerCircle;++j)
142 {
143 ScalarType X = cos(j*HorizontalAngle)*HorizRadius;
144 ScalarType Y = sin(j*HorizontalAngle)*HorizRadius;
145 NN.push_back(Point3<ScalarType>(X,Y,Z));
146 }
147 }
148 NN.push_back(Point3<ScalarType>(0,0,-1.0));
149 }
150
RecursiveOctahedron(int vn,std::vector<Point3<ScalarType>> & NN)151 static void RecursiveOctahedron(int vn, std::vector<Point3<ScalarType > > &NN)
152 {
153 OctaLevel pp;
154
155 int ll=10;
156 while(pow(4.0f,ll)+2>vn) ll--;
157
158 pp.Init(ll);
159 sort(pp.v.begin(),pp.v.end());
160 int newsize = unique(pp.v.begin(),pp.v.end())-pp.v.begin();
161 pp.v.resize(newsize);
162
163 NN=pp.v;
164 //Perturb(NN);
165 }
166
Perturb(std::vector<Point3<ScalarType>> & NN)167 static void Perturb(std::vector<Point3<ScalarType > > &NN)
168 {
169 float width=0.2f/sqrt(float(NN.size()));
170
171 typename std::vector<Point3<ScalarType> >::iterator vi;
172 for(vi=NN.begin(); vi!=NN.end();++vi)
173 {
174 Point3x pp(((float)rand())/RAND_MAX,
175 ((float)rand())/RAND_MAX,
176 ((float)rand())/RAND_MAX);
177 pp=pp*2.0-Point3x(1,1,1);
178 pp*=width;
179 (*vi)+=pp;
180 (*vi).Normalize();
181 }
182
183 }
184
185 /*
186 Trova la normale piu vicina a quella data.
187 Assume che tutte normale in ingresso sia normalizzata;
188 */
BestMatchingNormal(const Point3x & n,std::vector<Point3x> & nv)189 static int BestMatchingNormal(const Point3x &n, std::vector<Point3x> &nv)
190 {
191 int ret=-1;
192 ScalarType bestang=-1;
193 ScalarType cosang;
194 typename std::vector<Point3x>::iterator ni;
195 for(ni=nv.begin();ni!=nv.end();++ni)
196 {
197 cosang=(*ni).dot(n);
198 if(cosang>bestang) {
199 bestang=cosang;
200 ret=ni-nv.begin();
201 }
202 }
203 assert(ret>=0 && ret <int(nv.size()));
204 return ret;
205 }
206
207
208 private :
209 class OctaLevel
210 {
211 public:
212 std::vector<Point3x> v;
213 int level;
214 int sz;
215 int sz2;
216
Val(int i,int j)217 Point3x &Val(int i, int j) {
218
219 assert(i>=-sz2 && i<=sz2);
220 assert(j>=-sz2 && j<=sz2);
221 return v[i+sz2 +(j+sz2)*sz];
222 }
223 /*
224 * Only the first quadrant is generated and replicated onto the other ones.
225 *
226 * o lev == 1
227 * | \ sz2 = 2^lev = 2
228 * o - o sz = 5 (eg. all the points lie in a 5x5 squre)
229 * | \ | \
230 * o - o - o
231 *
232 * |
233 * V
234 *
235 * o
236 * | \ lev == 1
237 * o - o sz2 = 4
238 * | \ | \ sz = 9 (eg. all the points lie in a 9x9 squre)
239 * o - o - o
240 * | \ | \ | \
241 * o - o - o - o
242 * | \ | \ | \ | \
243 * o - o - o - o - o
244 *
245 *
246 */
Init(int lev)247 void Init(int lev)
248 {
249 sz2=pow(2.0f,lev);
250 sz=sz2*2+1;
251 v.resize(sz*sz,Point3x(0,0,0));
252 if(lev==0)
253 {
254 Val( 0,0)=Point3x( 0, 0, 1);
255 Val( 1,0)=Point3x( 1, 0, 0);
256 Val( 0,1)=Point3x( 0, 1, 0);
257 }
258 else
259 {
260 OctaLevel tmp;
261 tmp.Init(lev-1);
262 int i,j;
263 for(i=0;i<=sz2;++i)
264 for(j=0;j<=(sz2-i);++j)
265 {
266 if((i%2)==0 && (j%2)==0)
267 Val(i,j)=tmp.Val(i/2,j/2);
268 if((i%2)!=0 && (j%2)==0)
269 Val(i,j)=(tmp.Val((i-1)/2,j/2)+tmp.Val((i+1)/2,j/2))/2.0;
270 if((i%2)==0 && (j%2)!=0)
271 Val(i,j)=(tmp.Val(i/2,(j-1)/2)+tmp.Val(i/2,(j+1)/2))/2.0;
272 if((i%2)!=0 && (j%2)!=0)
273 Val(i,j)=(tmp.Val((i-1)/2,(j+1)/2)+tmp.Val((i+1)/2,(j-1)/2))/2.0;
274
275 Val( sz2-j, sz2-i)[0] = Val(i,j)[0]; Val( sz2-j, sz2-i)[1] = Val(i,j)[1]; Val( sz2-j, sz2-i)[2] = -Val(i,j)[2];
276 Val(-sz2+j, sz2-i)[0] =-Val(i,j)[0]; Val(-sz2+j, sz2-i)[1] = Val(i,j)[1]; Val(-sz2+j, sz2-i)[2] = -Val(i,j)[2];
277 Val( sz2-j,-sz2+i)[0] = Val(i,j)[0]; Val( sz2-j,-sz2+i)[1] =-Val(i,j)[1]; Val( sz2-j,-sz2+i)[2] = -Val(i,j)[2];
278 Val(-sz2+j,-sz2+i)[0] =-Val(i,j)[0]; Val(-sz2+j,-sz2+i)[1] =-Val(i,j)[1]; Val(-sz2+j,-sz2+i)[2] = -Val(i,j)[2];
279
280 Val(-i,-j)[0] = -Val(i,j)[0]; Val(-i,-j)[1] = -Val(i,j)[1]; Val(-i,-j)[2] = Val(i,j)[2];
281 Val( i,-j)[0] = Val(i,j)[0]; Val( i,-j)[1] = -Val(i,j)[1]; Val( i,-j)[2] = Val(i,j)[2];
282 Val(-i, j)[0] = -Val(i,j)[0]; Val(-i, j)[1] = Val(i,j)[1]; Val(-i, j)[2] = Val(i,j)[2];
283 }
284
285 typename std::vector<Point3<ScalarType> >::iterator vi;
286 for(vi=v.begin(); vi!=v.end();++vi)
287 (*vi).Normalize();
288 }
289 }
290 };
291 };
292 }
293 #endif
294