manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Added Fibonacci sampling, renamed to more meaningful names the sampling algs

This commit is contained in:
Paolo Cignoni 2015-11-04 23:49:35 +00:00
parent 1480d19996
commit 0f05ee423d
1 changed files with 84 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

110
vcg/math/gen_normal.h
View File

@ -55,6 +55,33 @@ static void Random(int vn, std::vector<Point3<ScalarType > > &NN)
} }
} }
} }
static Point3x FibonacciPt(int i, int n)
{
const ScalarType Phi = ScalarType(sqrt(5)*0.5 + 0.5);
const ScalarType phi = 2.0*M_PI* (i/Phi - floor(i/Phi));
ScalarType cosTheta = 1.0 - (2*i + 1.0)/ScalarType(n);
float sinTheta = 1 - cosTheta*cosTheta;
sinTheta = sqrt(std::min(ScalarType(1),std::max(ScalarType(0),sinTheta)));
return Point3x(
cos(phi)*sinTheta,
sin(phi)*sinTheta,
cosTheta);
}
// Implementation of the Spherical Fibonacci Point Sets
// according to the description of
// Spherical Fibonacci Mapping
// Benjamin Keinert, Matthias Innmann, Michael Sanger, Marc Stamminger
// TOG 2015
static void Fibonacci(int n, std::vector<Point3x > &NN)
{
NN.resize(n);
for(int i=0;i<n;++i)
NN[i]=FibonacciPt(i,n);
}
static void UniformCone(int vn, std::vector<Point3<ScalarType > > &NN, ScalarType AngleRad, Point3x dir=Point3x(0,1,0)) static void UniformCone(int vn, std::vector<Point3<ScalarType > > &NN, ScalarType AngleRad, Point3x dir=Point3x(0,1,0))
{ {
std::vector<Point3<ScalarType > > NNT; std::vector<Point3<ScalarType > > NNT;
@ -66,14 +93,14 @@ static void UniformCone(int vn, std::vector<Point3<ScalarType > > &NN, ScalarTyp
ScalarType Ratio = CapArea / (4.0*M_PI ); ScalarType Ratio = CapArea / (4.0*M_PI );
printf("----------AngleRad %f Angledeg %f ratio %f vn %i vn2 %i \n",AngleRad,math::ToDeg(AngleRad),Ratio,vn,int(vn/Ratio)); printf("----------AngleRad %f Angledeg %f ratio %f vn %i vn2 %i \n",AngleRad,math::ToDeg(AngleRad),Ratio,vn,int(vn/Ratio));
Uniform(vn/Ratio,NNT); Fibonacci(vn/Ratio,NNT);
printf("asked %i got %i (expecting %i instead of %i)\n", int(vn/Ratio), NNT.size(), int(NNT.size()*Ratio), vn); printf("asked %i got %i (expecting %i instead of %i)\n", int(vn/Ratio), NNT.size(), int(NNT.size()*Ratio), vn);
typename std::vector<Point3<ScalarType> >::iterator vi; typename std::vector<Point3<ScalarType> >::iterator vi;
ScalarType DotProd = cos(AngleRad); ScalarType cosAngle = cos(AngleRad);
for(vi=NNT.begin();vi!=NNT.end();++vi) for(vi=NNT.begin();vi!=NNT.end();++vi)
{ {
if(dir.dot(*vi) >= DotProd) NN.push_back(*vi); if(dir.dot(*vi) >= cosAngle) NN.push_back(*vi);
} }
} }
@ -90,7 +117,7 @@ static void UniformCone(int vn, std::vector<Point3<ScalarType > > &NN, ScalarTyp
// 2+ 2N*Sum(i=1 to N-1) sin(pi i/N). // 2+ 2N*Sum(i=1 to N-1) sin(pi i/N).
// The closed form of this summation // The closed form of this summation
// 2.0 - ( (2.0*N * sin (M_PI/N))/(cos(M_PI/N) - 1.0)); // 2.0 - ( (2.0*N * sin (M_PI/N))/(cos(M_PI/N) - 1.0));
static void Regular(int vn, std::vector<Point3<ScalarType > > &NN) static void DiscoBall(int vn, std::vector<Point3<ScalarType > > &NN)
{ {
// Guess the right N // Guess the right N
ScalarType N=0; ScalarType N=0;
@ -98,11 +125,9 @@ static void Regular(int vn, std::vector<Point3<ScalarType > > &NN)
for(N=1;N<vn;++N) for(N=1;N<vn;++N)
{ {
ScalarType expectedPoints = 2.0 - ( (2.0*N * sin (M_PI/N))/(cos(M_PI/N) - 1.0)); ScalarType expectedPoints = 2.0 - ( (2.0*N * sin (M_PI/N))/(cos(M_PI/N) - 1.0));
qDebug("N %f -> %f",N,expectedPoints);
if(expectedPoints >= vn) break; if(expectedPoints >= vn) break;
} }
ScalarType VerticalAngle = M_PI / N; ScalarType VerticalAngle = M_PI / N;
NN.push_back(Point3<ScalarType>(0,0,1.0)); NN.push_back(Point3<ScalarType>(0,0,1.0));
for (int i =1; i<N; ++i) for (int i =1; i<N; ++i)
@ -123,7 +148,7 @@ static void Regular(int vn, std::vector<Point3<ScalarType > > &NN)
NN.push_back(Point3<ScalarType>(0,0,-1.0)); NN.push_back(Point3<ScalarType>(0,0,-1.0));
} }
static void Uniform(int vn, std::vector<Point3<ScalarType > > &NN) static void RecursiveOctahedron(int vn, std::vector<Point3<ScalarType > > &NN)
{ {
OctaLevel pp; OctaLevel pp;
@ -136,7 +161,7 @@ static void Uniform(int vn, std::vector<Point3<ScalarType > > &NN)
pp.v.resize(newsize); pp.v.resize(newsize);
NN=pp.v; NN=pp.v;
Perturb(NN); //Perturb(NN);
} }
static void Perturb(std::vector<Point3<ScalarType > > &NN) static void Perturb(std::vector<Point3<ScalarType > > &NN)
@ -187,46 +212,79 @@ class OctaLevel
std::vector<Point3x> v; std::vector<Point3x> v;
int level; int level;
int sz; int sz;
int sz2;
Point3x &Val(int i, int j) { Point3x &Val(int i, int j) {
assert(i>=0 && i<sz);
assert(j>=0 && j<sz);
return v[i+j*sz];
}
assert(i>=-sz2 && i<=sz2);
assert(j>=-sz2 && j<=sz2);
return v[i+sz2 +(j+sz2)*sz];
}
/*
* Only the first quadrant is generated and replicated onto the other ones.
*
* o lev == 1
* | \ sz2 = 2^lev = 2
* o - o sz = 5 (eg. all the points lie in a 5x5 squre)
* | \ | \
* o - o - o
*
* |
* V
*
* o
* | \ lev == 1
* o - o sz2 = 4
* | \ | \ sz = 9 (eg. all the points lie in a 9x9 squre)
* o - o - o
* | \ | \ | \
* o - o - o - o
* | \ | \ | \ | \
* o - o - o - o - o
*
*
*/
void Init(int lev) void Init(int lev)
{ {
sz=pow(2.0f,lev+1)+1; sz2=pow(2.0f,lev);
v.resize(sz*sz); sz=sz2*2+1;
v.resize(sz*sz,Point3x(0,0,0));
if(lev==0) if(lev==0)
{ {
Val(0,0)=Point3x( 0, 0,-1); Val(0,1)=Point3x( 0, 1, 0); Val(0,2)=Point3x( 0, 0,-1); Val( 0,0)=Point3x( 0, 0, 1);
Val( 1,0)=Point3x( 1, 0, 0);
Val(1,0)=Point3x(-1, 0, 0); Val(1,1)=Point3x( 0, 0, 1); Val(1,2)=Point3x( 1, 0, 0); Val( 0,1)=Point3x( 0, 1, 0);
Val(2,0)=Point3x( 0, 0,-1); Val(2,1)=Point3x( 0,-1, 0); Val(2,2)=Point3x( 0, 0,-1);
} }
else else
{ {
OctaLevel tmp; OctaLevel tmp;
tmp.Init(lev-1); tmp.Init(lev-1);
int i,j; int i,j;
for(i=0;i<sz;++i) for(i=0;i<=sz2;++i)
for(j=0;j<sz;++j) for(j=0;j<=(sz2-i);++j)
{ {
if((i%2)==0 && (j%2)==0) if((i%2)==0 && (j%2)==0)
Val(i,j)=tmp.Val(i/2,j/2); Val(i,j)=tmp.Val(i/2,j/2);
if((i%2)!=0 && (j%2)==0) if((i%2)!=0 && (j%2)==0)
Val(i,j)=(tmp.Val(i/2+0,j/2)+tmp.Val(i/2+1,j/2))/2.0; Val(i,j)=(tmp.Val((i-1)/2,j/2)+tmp.Val((i+1)/2,j/2))/2.0;
if((i%2)==0 && (j%2)!=0) if((i%2)==0 && (j%2)!=0)
Val(i,j)=(tmp.Val(i/2,j/2+0)+tmp.Val(i/2,j/2+1))/2.0; Val(i,j)=(tmp.Val(i/2,(j-1)/2)+tmp.Val(i/2,(j+1)/2))/2.0;
if((i%2)!=0 && (j%2)!=0) if((i%2)!=0 && (j%2)!=0)
Val(i,j)=(tmp.Val(i/2+0,j/2+0)+tmp.Val(i/2+0,j/2+1)+tmp.Val(i/2+1,j/2+0)+tmp.Val(i/2+1,j/2+1))/4.0; Val(i,j)=(tmp.Val((i-1)/2,(j+1)/2)+tmp.Val((i+1)/2,(j-1)/2))/2.0;
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];
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];
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];
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];
Val(-i,-j)[0] = -Val(i,j)[0]; Val(-i,-j)[1] = -Val(i,j)[1]; Val(-i,-j)[2] = Val(i,j)[2];
Val( i,-j)[0] = Val(i,j)[0]; Val( i,-j)[1] = -Val(i,j)[1]; Val( i,-j)[2] = Val(i,j)[2];
Val(-i, j)[0] = -Val(i,j)[0]; Val(-i, j)[1] = Val(i,j)[1]; Val(-i, j)[2] = Val(i,j)[2];
} }
typename std::vector<Point3<ScalarType> >::iterator vi;
typename std::vector<Point3<ScalarType> >::iterator vi;
for(vi=v.begin(); vi!=v.end();++vi) for(vi=v.begin(); vi!=v.end();++vi)
(*vi).Normalize(); (*vi).Normalize();
} }
} }
}; };