More changes to the interface

This commit is contained in:
Nico Kruithof 2006-05-03 13:28:33 +00:00
parent 94c2aab68d
commit f0bf17c1a0
4 changed files with 240 additions and 243 deletions

View File

@ -1,5 +1,3 @@
// NGHK: remove later
#define CGAL_PROFILE
// examples/Skin_surface_3/skin_surface_simple.C
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Skin_surface_3.h>

View File

@ -24,13 +24,6 @@
CGAL_BEGIN_NAMESPACE
/////////////////////////////////////////////////
// //
// Virtual base class of a quadratic surface //
// This implementation is specific for //
// skin surfaces //
// //
/////////////////////////////////////////////////
template < class Kernel >
class Skin_surface_quadratic_surface_3 {
public:
@ -41,9 +34,32 @@ public:
typedef typename K::RT RT;
typedef Weighted_point<Point, RT> Weighted_point;
Skin_surface_quadratic_surface_3() {
Skin_surface_quadratic_surface_3(RT Qinput[], Point p, RT c)
: p(p), c(c)
{
for (int i=0; i<6; i++) Q[i] = Qinput[i];
}
virtual ~Skin_surface_quadratic_surface_3() {};
template <class Input_point>
RT value(Input_point const &x) const {
typedef Cartesian_converter<typename Input_point::R, K> Converter;
Vector v = Converter()(x) - p;
return
v.x()*(Q[0]*v.x()) +
v.y()*(Q[1]*v.x()+Q[2]*v.y()) +
v.z()*(Q[3]*v.x()+Q[4]*v.y()+Q[5]*v.z()) +
c;
}
Vector gradient(Point const &x) {
std::cout << "NGHK: NOT YET IMPLEMENTED" << std::endl;
// NGHK: TODO:
return (x-p);
}
/// Construct the intersection point with the segment (p0,p1)
template <class Input_point>
Input_point to_surface(const Input_point &p0, const Input_point &p1) {
@ -69,247 +85,215 @@ public:
}
return mid;
};
// virtual Point to_surface(Point const &p, Vector const &v) = 0;
inline Point to_surface(const Segment &s) {
return to_surface(s.source(), s.target());
}
// // Gradient descent
// virtual Point to_surface(Point const &p0) = 0;
// compute the function value of p
template <class Input_point>
RT value(const Input_point &p) const {
typedef Cartesian_converter<typename Input_point::R, K> Converter;
return _value(Converter()(p));
}
virtual RT _value(const Point &p) const = 0;
// Compute the gradient in p
virtual Vector gradient(const Point &p) = 0;
// return the dimension of the delaunay simplex:
virtual int dimension() const = 0;
// return a continuous density function on the skin surface:
virtual RT sq_density(const Point &p) = 0;
};
/////////////////////////////////////////////////
// //
// Sphere //
// //
// orient*|x-p|^2/s - P == 0 //
// orient*|x-p|^2/s - P == 0 //
// //
// p: center of the sphere //
// P: squared radius of the sphere //
// orient: orientation of the sphere //
// //
/////////////////////////////////////////////////
template < class Kernel >
class Skin_surface_sphere_3 : public Skin_surface_quadratic_surface_3<Kernel> {
public:
typedef Skin_surface_quadratic_surface_3<Kernel> Parent;
typedef typename Parent::Point Point;
typedef typename Parent::Vector Vector;
typedef typename Parent::RT RT;
typedef typename Parent::Weighted_point Weighted_point;
Skin_surface_sphere_3(Weighted_point wp, RT s, int orient)
: Skin_surface_quadratic_surface_3<Kernel>(),
wp(wp), s(s), orient(orient) {
assert((orient == -1) || (orient == 1));
}
RT _value(Point const &x) const {
return orient * (squared_distance(x, wp)/s - wp.weight());
}
Point to_surface(const Point &p, const Vector &v) {
Vector pc = p-wp;
RT sq_d = v*v;
RT top = -pc*v/sq_d;
RT extr_val = squared_distance(p+top*v, wp) - s*wp.weight();
if (extr_val > 0) {
// std::cerr << "im. intersection[" << dimension() << "] "
// << extr_val << "\n";
return p;
}
RT d = sqrt(-extr_val/sq_d);
// t should be in [0,1]
RT t, t1;
t = top + d; t1 = top - d;
if (t1*t1 < t*t) t = t1;
CGAL_assertion(abs(value(p + t*v)) < 0.001);
return p + t*v;
}
Point to_surface(Point const &p) {
Vector pc = p-wp;
return wp + sqrt(s*wp.weight()/(pc*pc))*pc;
}
Vector gradient(Point const &p) {
return orient*(p-wp);
}
int dimension() const {
if (orient == 1) return 0; else return 3;
}
RT sq_density(Point const &p) {
assert(wp.weight() > 0);
return s*wp.weight();
}
private:
Weighted_point wp;
RT s;
int orient;
RT Q[6];
Point p;
RT c;
};
/* *************************************************
* Hyperboloid
*
* Point p.
* Let y = (p-wp)*t/|t|
* x = |p-wp|^2 - y^2
*
* orient*((1-s) x^2 - s y^2 + s(1-s) W = 0
* orient*((1-s)|p-wp|^2 - y^2 + s(1-s)W = 0
*
***************************************************/
// template < class Kernel >
// class Skin_surface_quadratic_surface_3 {
// public:
// typedef Kernel K;
// typedef typename K::Point_3 Point;
// typedef typename K::Vector_3 Vector;
// typedef typename K::Segment_3 Segment;
// typedef typename K::RT RT;
// typedef Weighted_point<Point, RT> Weighted_point;
template < class Kernel >
class Skin_surface_hyperboloid_3 : public Skin_surface_quadratic_surface_3<Kernel> {
public:
typedef Skin_surface_quadratic_surface_3<Kernel> Parent;
typedef typename Parent::Point Point;
typedef typename Parent::Vector Vector;
typedef typename Parent::RT RT;
typedef typename Parent::Weighted_point Weighted_point;
// Skin_surface_quadratic_surface_3() {
// }
// virtual ~Skin_surface_quadratic_surface_3() {};
// /// Construct the intersection point with the segment (p0,p1)
// template <class Input_point>
// Input_point to_surface(const Input_point &p0, const Input_point &p1) {
// // NGHK: We use the kernel trick again: DOCUMENT!!!!
// typedef Cartesian_converter<typename Input_point::R, K> Converter;
// Converter conv;
Skin_surface_hyperboloid_3(Weighted_point wp, Vector t, RT s, int orient)
: Skin_surface_quadratic_surface_3<Kernel>(), wp(wp), t(t), s(s), orient(orient) {
assert((orient == -1) || (orient == 1));
sq_t = t*t;
}
// Input_point pp0 = p0, pp1 = p1, mid;
// double sq_d = to_double(squared_distance(pp0,pp1));
// if (value(conv(pp1)) < value(conv(pp0))) {
// std::swap(pp0, pp1);
// }
// while (sq_d > 1.e-10) {
// mid = midpoint(pp0,pp1);
// if (value(conv(mid)) > 0) {
// pp1 = mid;
// } else {
// pp0 = mid;
// }
// sq_d /= 4;
// }
// return mid;
// };
// // virtual Point to_surface(Point const &p, Vector const &v) = 0;
// inline Point to_surface(const Segment &s) {
// return to_surface(s.source(), s.target());
// }
// // // Gradient descent
// // virtual Point to_surface(Point const &p0) = 0;
// // compute the function value of p
// template <class Input_point>
// RT value(const Input_point &p) const {
// typedef Cartesian_converter<typename Input_point::R, K> Converter;
// return _value(Converter()(p));
// }
// virtual RT _value(const Point &p) const = 0;
// // Compute the gradient in p
// virtual Vector gradient(const Point &p) = 0;
// // // NGHK: REMOVE:
// virtual int dimension() const = 0;
// // // return a continuous density function on the skin surface:
// // virtual RT sq_density(const Point &p) = 0;
// };
// /* *************************************************
// * Hyperboloid
// *
// * Point p.
// * Let y = (p-wp)*t/|t|
// * x = |p-wp|^2 - y^2
// *
// * orient*((1-s) x^2 - s y^2 + s(1-s) W = 0
// * orient*((1-s)|p-wp|^2 - y^2 + s(1-s)W = 0
// *
// ***************************************************/
// template < class Kernel >
// class Skin_surface_hyperboloid_3 : public Skin_surface_quadratic_surface_3<Kernel> {
// public:
// typedef Skin_surface_quadratic_surface_3<Kernel> Parent;
// typedef typename Parent::Point Point;
// typedef typename Parent::Vector Vector;
// typedef typename Parent::RT RT;
// typedef typename Parent::Weighted_point Weighted_point;
// Skin_surface_hyperboloid_3(Weighted_point wp, Vector t, RT s, int orient)
// : Skin_surface_quadratic_surface_3<Kernel>(), wp(wp), t(t), s(s), orient(orient) {
// assert((orient == -1) || (orient == 1));
// sq_t = t*t;
// }
RT _value(Point const &x) const {
Vector dir = x-wp;
RT tmp = dir*t;
tmp = tmp*tmp/sq_t;
// RT _value(Point const &x) const {
// Vector dir = x-wp;
// RT tmp = dir*t;
// tmp = tmp*tmp/sq_t;
return orient * (dir*dir/s - tmp/(s*(1-s))) + wp.weight();
}
// Point to_surface(Point const &p0, Point const &p1) {
// assert(value(p0) * value(p1) <= 0);
// return orient * (dir*dir/s - tmp/(s*(1-s))) + wp.weight();
// }
// // Point to_surface(Point const &p0, Point const &p1) {
// // assert(value(p0) * value(p1) <= 0);
// Vector p0c = p0-wp;
// Vector p1p0 = p1-p0;
// RT sq_d = p1p0*p1p0;
// // Vector p0c = p0-wp;
// // Vector p1p0 = p1-p0;
// // RT sq_d = p1p0*p1p0;
// RT p0_sym = p0c*t;
// // RT p0_sym = p0c*t;
// // RT p1_sym = (p1-wp)*t;
// // RT d_sym = p0_sym-p1_sym;
// // RT den = ((1-s)*sq_d - d_sym*d_sym/sq_t);
// // RT top = -((1-s)*(p0c*p1p0) + p0_sym*d_sym/sq_t) / den;
// // RT extr_val =
// // orient*((1-s)*p0c*p0c-p0_sym*p0_sym/sq_t-s*(1-s)*wp.weight());
// // {
// // Point extr = p0+top*p1p0;
// // Vector dir = extr-wp;
// // RT tmp = dir*t; tmp *= tmp/sq_t;
// // extr_val = orient*(-tmp + (1-s)*dir*dir + s*(1-s)*wp.weight());
// // }
// // RT d = sqrt(-orient*extr_val/den);
// // RT t, t1;
// // t = top + d; t1 = top - d;
// // if ((2*t1-1)*(2*t1-1) < (2*t-1)*(2*t-1)) t = t1;
// // if (t < 0) {
// // std::cerr << "Hl[" << dimension() <<"] " << t << "\n";
// // return p0;
// // } else if (t > 1) {
// // std::cerr << "Hh[" << dimension() <<"] " << t << "\n";
// // return p1;
// // } else {
// // assert (std::abs(value(p0 + t*p1p0)) < 0.001);
// // return p0 + t*p1p0;
// // }
// // }
// Point to_surface(Point const &p, Vector const &v){
// Vector pc = p-wp;
// Point p1 = p + v;
// RT sq_d = v*v;
// RT p_sym = pc*t;
// RT p1_sym = (p1-wp)*t;
// RT d_sym = p0_sym-p1_sym;
// RT d_sym = p_sym-p1_sym;
// RT den = ((1-s)*sq_d - d_sym*d_sym/sq_t);
// RT top = -((1-s)*(p0c*p1p0) + p0_sym*d_sym/sq_t) / den;
// RT top = -((1-s)*(pc*v) + p_sym*d_sym/sq_t) / den;
// RT extr_val =
// orient*((1-s)*p0c*p0c-p0_sym*p0_sym/sq_t-s*(1-s)*wp.weight());
// orient*((1-s)*pc*pc-p_sym*p_sym/sq_t-s*(1-s)*wp.weight());
// {
// Point extr = p0+top*p1p0;
// Point extr = p+top*v;
// Vector dir = extr-wp;
// RT tmp = dir*t; tmp *= tmp/sq_t;
// extr_val = orient*(-tmp + (1-s)*dir*dir + s*(1-s)*wp.weight());
// }
// RT d = sqrt(-orient*extr_val/den);
// RT d = -orient*extr_val/den;
// if (d<0) return p;
// d = sqrt(d);
// RT t, t1;
// t = top + d; t1 = top - d;
// if ((2*t1-1)*(2*t1-1) < (2*t-1)*(2*t-1)) t = t1;
// if (t1*t1 < t*t) t = t1;
// if (t < 0) {
// std::cerr << "Hl[" << dimension() <<"] " << t << "\n";
// return p0;
// } else if (t > 1) {
// std::cerr << "Hh[" << dimension() <<"] " << t << "\n";
// return p1;
// } else {
// assert (std::abs(value(p0 + t*p1p0)) < 0.001);
// return p0 + t*p1p0;
// }
// assert (std::abs(value(p + t*v)) < 0.001);
// return p + t*v;
// }
Point to_surface(Point const &p, Vector const &v){
Vector pc = p-wp;
Point p1 = p + v;
RT sq_d = v*v;
// Point to_surface(Point const &p0) {
// return to_surface(p0, gradient(p0));
// }
// Vector gradient(Point const &p) {
// // -s x + (1-s) y
// Vector v = p - wp;
// Vector vt = (v*t)/(t*t)*t;
// return orient*((1-s)*v - vt);
// }
// int dimension() const {
// if (orient == 1) return 1; else return 2;
// }
// RT sq_density(Point const &p) {
// Vector v = p - wp;
// RT vt = v*t;
// RT scale = 1-s;
// if (s < RT(.5)) {
// scale = (scale-s)/(scale*scale);
// } else {
// scale = (s-scale)/(s*s);
// }
// assert(scale >= 0);
RT p_sym = pc*t;
RT p1_sym = (p1-wp)*t;
RT d_sym = p_sym-p1_sym;
RT den = ((1-s)*sq_d - d_sym*d_sym/sq_t);
RT top = -((1-s)*(pc*v) + p_sym*d_sym/sq_t) / den;
RT extr_val =
orient*((1-s)*pc*pc-p_sym*p_sym/sq_t-s*(1-s)*wp.weight());
{
Point extr = p+top*v;
Vector dir = extr-wp;
RT tmp = dir*t; tmp *= tmp/sq_t;
extr_val = orient*(-tmp + (1-s)*dir*dir + s*(1-s)*wp.weight());
}
RT d = -orient*extr_val/den;
if (d<0) return p;
d = sqrt(d);
RT t, t1;
t = top + d; t1 = top - d;
if (t1*t1 < t*t) t = t1;
assert (std::abs(value(p + t*v)) < 0.001);
return p + t*v;
}
Point to_surface(Point const &p0) {
return to_surface(p0, gradient(p0));
}
Vector gradient(Point const &p) {
// -s x + (1-s) y
Vector v = p - wp;
Vector vt = (v*t)/(t*t)*t;
return orient*((1-s)*v - vt);
}
int dimension() const {
if (orient == 1) return 1; else return 2;
}
RT sq_density(Point const &p) {
Vector v = p - wp;
RT vt = v*t;
RT scale = 1-s;
if (s < RT(.5)) {
scale = (scale-s)/(scale*scale);
} else {
scale = (s-scale)/(s*s);
}
assert(scale >= 0);
return squared_distance(wp, p) - scale*vt*vt/(t*t);
}
private:
Weighted_point wp;
Vector t;
RT s, sq_t;
int orient;
// return squared_distance(wp, p) - scale*vt*vt/(t*t);
// }
// private:
// Weighted_point wp;
// Vector t;
// RT s, sq_t;
// int orient;
};
// };
CGAL_END_NAMESPACE

View File

@ -67,9 +67,9 @@ public:
// typedef typename Triangulated_mixed_complex_traits::RT TMC_RT;
// typedef Weighted_point<TMC_Point,TMC_RT> TMC_Weighted_point;
// typedef Skin_surface_quadratic_surface_3<Polyhedron_traits> QuadrSurface;
typedef Skin_surface_sphere_3<Surface_traits> Sphere_surface;
typedef Skin_surface_hyperboloid_3<Surface_traits> Hyperboloid_surface;
typedef Skin_surface_quadratic_surface_3<Surface_traits> Quadratic_surface;
// typedef Skin_surface_sphere_3<Surface_traits> Sphere_surface;
// typedef Skin_surface_hyperboloid_3<Surface_traits> Hyperboloid_surface;
typedef Weighted_converter_3<
Cartesian_converter < typename Regular_traits::Bare_point::R,
@ -111,8 +111,9 @@ public:
{
vh = s;
create_sphere(r2s_converter(vh->point().point()),
r2s_converter(vh->point().weight()),
shrink);
-r2s_converter(vh->point().weight()),
shrink,
1);
break;
}
case 1:
@ -170,7 +171,8 @@ public:
r2s_converter(ch->vertex(1)->point()),
r2s_converter(ch->vertex(2)->point()),
r2s_converter(ch->vertex(3)->point())),
shrink-1);
1-shrink,
-1);
}
}
}
@ -183,22 +185,34 @@ public:
void create_sphere(const Surface_point &c,
const Surface_RT &w,
const Surface_RT &s) {
if (s >= 0) {
surf = new Sphere_surface(Surface_weighted_point(c,w), s, 1);
} else {
surf = new Sphere_surface(Surface_weighted_point(c,w), -s, -1);
}
const Surface_RT &s,
const int orient) {
Q[1] = Q[3] = Q[4] = 0;
Q[0] = Q[2] = Q[5] = orient*(1-s);
surf = new Quadratic_surface(Q, c, s*(1-s)*w);
}
void create_hyperboloid(const Surface_point &c,
const Surface_RT &w,
const Surface_vector &v,
const Surface_vector &t,
const Surface_RT &s,
int orient) {
surf = new Hyperboloid_surface(Surface_weighted_point(c,w), v, s, orient);
const int orient) {
Surface_RT den = t*t;
Q[0] = orient*(- t.x()*t.x()/den + (1-s));
Q[1] = orient*(-2*t.y()*t.x()/den);
Q[2] = orient*(- t.y()*t.y()/den + (1-s));
Q[3] = orient*(-2*t.z()*t.x()/den);
Q[4] = orient*(-2*t.z()*t.y()/den);
Q[5] = orient*(- t.z()*t.z()/den + (1-s));
surf = new Quadratic_surface(Q, c, s*(1-s)*w);
}
Surface_RT Q[6];
R2S_converter r2s_converter;
};

View File

@ -53,7 +53,7 @@ typedef Triangulated_mixed_complex::Finite_vertices_iterator
#include <fstream>
int main(int argc, char *argv[]) {
double shrink = .85;
double shrink = .5;
std::list<Weighted_point> l;
while (argc>1) {
@ -71,6 +71,7 @@ int main(int argc, char *argv[]) {
for (Tmc_Finite_vertices_iterator vit = tmc.finite_vertices_begin();
vit != tmc.finite_vertices_end(); vit++) {
std::cout << std::endl;
if (tmc.is_infinite(vit->cell())) {
std::cerr << "ERROR: is_infinite (main)" << std::endl;
}
@ -82,7 +83,7 @@ int main(int argc, char *argv[]) {
cell != cells.end(); cell++) {
if (!tmc.is_infinite(*cell)) {
Quadratic_surface::RT val2 = (*cell)->surf->value(vit->point());
CGAL_assertion(val == val2);
CGAL_assertion(val == val2);
}
}
}