mirror of https://github.com/CGAL/cgal
324 lines
6.9 KiB
C++
324 lines
6.9 KiB
C++
// Copyright (c) 2003-2006 INRIA Sophia-Antipolis (France).
|
|
// All rights reserved.
|
|
//
|
|
// This file is part of CGAL (www.cgal.org).
|
|
// You can redistribute it and/or modify it under the terms of the GNU
|
|
// General Public License as published by the Free Software Foundation,
|
|
// either version 3 of the License, or (at your option) any later version.
|
|
//
|
|
// Licensees holding a valid commercial license may use this file in
|
|
// accordance with the commercial license agreement provided with the software.
|
|
//
|
|
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
|
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
|
//
|
|
// $URL$
|
|
// $Id$
|
|
//
|
|
//
|
|
// Author(s) : Steve OUDOT, Laurent Rineau
|
|
|
|
///////////////// Code for functions of famous surfaces /////////////////
|
|
|
|
|
|
// Sphere (r=1)
|
|
|
|
double sphere_function (double x, double y, double z) {
|
|
double x2=x*x, y2=y*y, z2=z*z;
|
|
|
|
return x2+y2+z2-1;
|
|
}
|
|
|
|
|
|
// Ellipsoid (r=1)
|
|
|
|
double ellipsoid_function (double x, double y, double z) {
|
|
double x2=x*x, y2=y*y, z2=z*z;
|
|
|
|
return x2+2*y2+4*z2-1;
|
|
}
|
|
|
|
|
|
// Torus (r=2)
|
|
|
|
double torus_function (double x, double y, double z) {
|
|
double x2=x*x, y2=y*y, z2=z*z;
|
|
double x4=x2*x2, y4=y2*y2, z4=z2*z2;
|
|
|
|
return x4 + y4 + z4 + 2 *x2* y2 + 2*
|
|
x2*z2 + 2*y2* z2 - 5 *x2 + 4* y2 - 5*z2+4;
|
|
}
|
|
|
|
|
|
// "Chair" (r=6)
|
|
|
|
double chair_function (double x, double y, double z) {
|
|
double x2=x*x, y2=y*y, z2=z*z;
|
|
double x4=x2*x2, y4=y2*y2, z4=z2*z2;
|
|
|
|
return x4-1.2*x2*y2+3.6*x2*z2-7.50*x2+y4+3.6*y2*z2-7.50*y2+.2*z4-7.50*z2+64.0625-16.0*z*y2+16.0*x2*z;
|
|
}
|
|
|
|
|
|
// "Tanglecube" (r=4)
|
|
|
|
double tanglecube_function (double x, double y, double z) {
|
|
double x2=x*x, y2=y*y, z2=z*z;
|
|
double x4=x2*x2, y4=y2*y2, z4=z2*z2;
|
|
|
|
return x4 - 5*x2 + y4 - 5*y2 + z4 - 5*z2 + 11.8;
|
|
}
|
|
|
|
double cube_function (double x, double y, double z){
|
|
if( x < 1 && -x < 1 &&
|
|
y < 1 && -y < 1 &&
|
|
z < 1 && -z < 1 )
|
|
return -1.;
|
|
return 1.;
|
|
}
|
|
|
|
// Barth's octic surface (degree 8)
|
|
|
|
double octic_function (double x, double y, double z) { // r=2
|
|
double x2=x*x, y2=y*y, z2=z*z;
|
|
double x4=x2*x2, y4=y2*y2, z4=z2*z2;
|
|
double x6=x4*x2, y6=y4*y2, z6=z4*z2;
|
|
double x8=x4*x4, y8=y4*y4, z8=z4*z4;
|
|
|
|
return 43.30495169 *x2* y2 + 43.30495169 *x2* z2 + 43.30495169 *y2 * z2 + 44.36067980 *x6* y2 + 44.36067980* x6 * z2 + 66.54101970* x4* y4 + 66.54101970* x4 * z4 + 44.36067980 *x2 * y6 - 11.70820393* x2 - 11.70820393* y2 - 11.70820393* z2 + 37.65247585 *x4 + 37.65247585 *y4 + 37.65247585* z4 + 11.09016995* x8 + 11.09016995 *y8 + 11.09016995* z8 + 133.0820394 *x2* y4* z2 + 133.0820394 *x2 * y2 * z4 + 44.36067980* x2 * z6 + 44.36067980 *y6 * z2 + 66.54101970 *y4 * z4 + 44.36067980 *y2 * z6 + 133.0820394* x4* y2 * z2 - 91.95742756 *x4 * y2 - 91.95742756 *x4 *z2 - 91.95742756* x2 * y4 - 91.95742756 *x2 *z4 - 91.95742756* y4 * z2 - 183.9148551* x2 *y2 *z2 - 30.65247585 *x6 - 30.65247585* y6 - 91.95742756 *y2 * z4 - 30.65247585* z6 + 1.618033988;
|
|
}
|
|
|
|
|
|
// "Heart"
|
|
|
|
double heart_function (double x, double y, double z) { // radius = 2
|
|
|
|
return (2*x*x+y*y+z*z-1)*(2*x*x+y*y+z*z-1)*(2*x*x+y*y+z*z-1) - (0.1*x*x+y*y)*z*z*z;
|
|
}
|
|
|
|
|
|
// Klein's bottle
|
|
|
|
double klein_function (double x, double y, double z) { // radius = 4
|
|
|
|
return (x*x+y*y+z*z+2*y-1)*((x*x+y*y+z*z-2*y-1)*(x*x+y*y+z*z-2*y-1)-8*z*z)+16*x*z*(x*x+y*y+z*z-2*y-1);
|
|
}
|
|
|
|
|
|
// Ring
|
|
|
|
double ring_function (double x, double y, double z) { // radius = ?
|
|
double e=0.1;
|
|
|
|
double f1 = x*x+y*y+z*z-1;
|
|
double f2 = x;
|
|
|
|
f1 = f1*f1-e*e;
|
|
f2 = f2*f2-e*e;
|
|
|
|
if (f1 < 0 && f2 < 0)
|
|
return -1.;
|
|
else if (f1 > 0 || f2 > 0)
|
|
return 1.;
|
|
else
|
|
return 0.;
|
|
}
|
|
|
|
|
|
// False knot
|
|
|
|
double false_knot_function (double x, double y, double z) { // radius = 1
|
|
double d=1.2, e=0.1;
|
|
|
|
double f1 = x*(x*x-z*z)-2*x*z*z-y*y+d*d-x*x-z*z-y*y;
|
|
|
|
double m2 = z*(x*x-z*z)+2*x*x*z;
|
|
double f2 = 4*y*y*(d*d-x*x-z*z-y*y) - m2*m2;
|
|
|
|
f1 = f1*f1-e*e;
|
|
f2 = f2*f2-e*e;
|
|
|
|
if (f1 < 0 && f2 < 0)
|
|
return -1.;
|
|
else if (f1 > 0 || f2 > 0)
|
|
return 1.;
|
|
else
|
|
return 0.;
|
|
}
|
|
|
|
|
|
// Knot 1
|
|
|
|
void puiss(double& x, double& y, int n) {
|
|
|
|
double xx = 1, yy = 0;
|
|
|
|
while(n>0) {
|
|
if (n&1) {
|
|
double xxx = xx, yyy = yy;
|
|
xx = xxx*x - yyy*y;
|
|
yy = xxx*y + yyy*x;
|
|
}
|
|
|
|
double xxx = x, yyy = y;
|
|
x=xxx*xxx-yyy*yyy;
|
|
y=2*xxx*yyy;
|
|
|
|
n/=2;
|
|
}
|
|
|
|
x = xx;
|
|
y = yy;
|
|
}
|
|
|
|
|
|
|
|
double knot1_function (double a, double b, double c) { // radius = 4
|
|
double e=0.1;
|
|
|
|
double x, y, z, t, den;
|
|
den=1+a*a+b*b+c*c;
|
|
x=2*a/den;
|
|
y=2*b/den;
|
|
z=2*c/den;
|
|
t=(1-a*a-b*b-c*c)/den;
|
|
|
|
double x3=x, y3=y, z2=z, t2=t;
|
|
puiss(x3,y3,3);
|
|
puiss(z2,t2,2);
|
|
|
|
double f1 = z2-x3;
|
|
double f2 = t2-y3;
|
|
|
|
f1 = f1*f1;
|
|
f2 = f2*f2;
|
|
e=e*e/(den-1);
|
|
|
|
if (f1+f2-e < 0)
|
|
return -1.;
|
|
else if (f1+f2-e > 0)
|
|
return 1.;
|
|
else
|
|
return 0.;
|
|
}
|
|
|
|
/*
|
|
double knot1_function (double x, double y, double z) { // radius = 4
|
|
double e=0.1;
|
|
|
|
double c1, c2, c3, c4, den;
|
|
den=1+x*x+y*y+z*z;
|
|
c1=2*x/den;
|
|
c2=2*y/den;
|
|
c3=2*z/den;
|
|
c4=(1-x*x-y*y-z*z)/den;
|
|
|
|
double f1 = c1*(c1*c1-c2*c2)-2*c1*c2*c2-c3*c3+c4*c4;
|
|
double f2 = c2*(c1*c1-c2*c2)+2*c1*c1*c2-2*c3*c4;
|
|
|
|
f1 = f1*f1;
|
|
f2 = f2*f2;
|
|
e=e*e/(den-1);
|
|
|
|
if (f1+f2-e < 0)
|
|
return -1.;
|
|
else if (f1+f2-e > 0)
|
|
return 1.;
|
|
else
|
|
return 0.;
|
|
}
|
|
*/
|
|
|
|
double knot1_surf1_function (double x, double y, double z) { // radius = 4
|
|
double c1, c2, c3, c4, den;
|
|
den=1+x*x+y*y+z*z;
|
|
c1=2*x/den;
|
|
c2=2*y/den;
|
|
c3=2*z/den;
|
|
c4=(1-x*x-y*y-z*z)/den;
|
|
|
|
return c1*(c1*c1-c2*c2)-2*c1*c2*c2-c3*c3+c4*c4;
|
|
}
|
|
|
|
|
|
double knot1_surf2_function (double x, double y, double z) { // radius = 4
|
|
double c1, c2, c3, c4, den;
|
|
den=1+x*x+y*y+z*z;
|
|
c1=2*x/den;
|
|
c2=2*y/den;
|
|
c3=2*z/den;
|
|
c4=(1-x*x-y*y-z*z)/den;
|
|
|
|
return c2*(c1*c1-c2*c2)+2*c1*c1*c2-2*c3*c4;
|
|
}
|
|
|
|
|
|
// Knot 2
|
|
|
|
double knot2_function (double a, double b, double c) { // radius = 4
|
|
double e=0.025;
|
|
|
|
double x, y, z, t, den;
|
|
den=1+a*a+b*b+c*c;
|
|
x=2*a/den;
|
|
y=2*b/den;
|
|
z=2*c/den;
|
|
t=(1-a*a-b*b-c*c)/den;
|
|
|
|
double x7=x, y7=y, x13=x, y13=y;
|
|
puiss(x7,y7,7);
|
|
puiss(x13,y13,13);
|
|
|
|
double z3=z, t3=t;
|
|
puiss(z3,t3,3);
|
|
|
|
double f1t = (z3-x7)*(z3-x7+100*x13) - (t3-y7)*(t3-y7+100*y13);
|
|
double f2t = (z3-x7)*(t3-y7+100*y13) + (t3-y7)*(z3-x7+100*x13);
|
|
double f1 = f1t*(z3-x7-100*x13) - f2t*(t3-y7-100*y13);
|
|
double f2 = f1t*(t3-y7-100*y13) + f2t*(z3-x7-100*x13);
|
|
|
|
f1 = f1*f1;
|
|
f2 = f2*f2;
|
|
e=e*e/(den-1);
|
|
|
|
if (f1+f2-e < 0)
|
|
return -1.;
|
|
else if (f1+f2-e > 0)
|
|
return 1.;
|
|
else
|
|
return 0.;
|
|
}
|
|
|
|
|
|
// Knot 3
|
|
|
|
|
|
double knot3_function (double a, double b, double c) { // radius = 4
|
|
double e=0.025;
|
|
|
|
double x, y, z, t, den;
|
|
den=1+a*a+b*b+c*c;
|
|
x=2*a/den;
|
|
y=2*b/den;
|
|
z=2*c/den;
|
|
t=(1-a*a-b*b-c*c)/den;
|
|
|
|
double x19=x, y19=y, z17=z, t17=t;
|
|
puiss(x19,y19,19);
|
|
puiss(z17,t17,17);
|
|
|
|
double f1 = z17-x19;
|
|
double f2 = t17-y19;
|
|
|
|
f1 = f1*f1;
|
|
f2 = f2*f2;
|
|
e=e*e/(den-1);
|
|
|
|
if (f1+f2-e < 0)
|
|
return -1.;
|
|
else if (f1+f2-e > 0)
|
|
return 1.;
|
|
else
|
|
return 0.;
|
|
}
|