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