--- a/.issues/c420501cc5285bbc Wed Mar 19 12:50:35 2008 -0400
+++ b/.issues/c420501cc5285bbc Wed Mar 19 13:50:51 2008 -0400
@@ -1,7 +1,7 @@
From artemis Tue Feb 26 10:25:07 2008
From: Dmitriy Morozov <morozov@cs.duke.edu>
Date: Tue, 26 Feb 2008 05:22:56 -0500
-State: open
+State: fixed
Subject: Non-optimized CGAL runtime error
Message-Id: <c420501cc5285bbc-0-artemis@metatron>
@@ -59,3 +59,13 @@
In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
state=open
+
+From artemis Wed Mar 19 17:17:47 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 13:17:30 -0400
+Subject: Fixed for ar-vineyard
+Message-Id: <c420501cc5285bbc-feb676745fc16b8c-artemis@metatron>
+References: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+
+Fixed for ar-vineyard.
--- a/examples/ar-vineyard/ar-simplex3d.h Wed Mar 19 12:50:35 2008 -0400
+++ b/examples/ar-vineyard/ar-simplex3d.h Wed Mar 19 13:50:51 2008 -0400
@@ -56,10 +56,11 @@
ARSimplex3D(const Edge& e);
ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices,
- Facet_circulator facet_bg);
+ const Delaunay& Dt, Facet_circulator facet_bg);
ARSimplex3D(const Facet& f);
- ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices);
+ ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices,
+ const Delaunay& Dt);
ARSimplex3D(const Cell& c, const Point& z);
--- a/examples/ar-vineyard/ar-simplex3d.hpp Wed Mar 19 12:50:35 2008 -0400
+++ b/examples/ar-vineyard/ar-simplex3d.hpp Wed Mar 19 13:50:51 2008 -0400
@@ -41,19 +41,15 @@
}
ARSimplex3D::
-ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, Facet_circulator facet_bg)
+ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, const Delaunay& Dt, Facet_circulator facet_bg)
{
Cell_handle c = e.first;
Parent::add(c->vertex(e.second));
Parent::add(c->vertex(e.third));
Facet_circulator cur = facet_bg;
+ while (Dt.is_infinite(*cur)) ++cur;
SimplexPhiMap::const_iterator cur_iter = simplices.find(ARSimplex3D(*cur));
- while (cur_iter == simplices.end())
- {
- ++cur;
- cur_iter = simplices.find(ARSimplex3D(*cur));
- }
RealValue min = cur_iter->first.alpha();
RealValue phi_const_min = cur_iter->first.phi_const();
@@ -69,10 +65,8 @@
int i0 = (*cur).first->index(*v1);
int i1 = (*cur).first->index(*v2);
int i = 6 - i0 - i1 - (*cur).second;
- Point p3 = (*cur).first->vertex(i)->point();
- cur_iter = simplices.find(ARSimplex3D(*cur));
- if (cur_iter == simplices.end()) // cur is infinite
+ if (Dt.is_infinite(cur->first->vertex(i)))
{
++cur; continue;
// FIXME: what do we do with infinite cofaces (i.e., what
@@ -80,9 +74,11 @@
// infinite?
}
+ Point p3 = (*cur).first->vertex(i)->point();
if (CGAL::side_of_bounded_sphere(p1, p2, p3) == CGAL::ON_BOUNDED_SIDE)
attached_ = true;
+ SimplexPhiMap::const_iterator cur_iter = simplices.find(ARSimplex3D(*cur));
RealValue val = cur_iter->first.alpha();
if (val < min) min = val;
RealValue phi_const_val = cur_iter->first.phi_const();
@@ -105,7 +101,7 @@
s_ = CGAL::squared_distance(z, K::Segment_3(p1,p2).supporting_line());
- Point origin;
+ Point origin(0,0,0);
Point cc = origin + ((p1 - origin) + (p2 - origin))/2; // CGAL is funny
v_ = CGAL::squared_distance(z, cc) - s_;
}
@@ -120,7 +116,7 @@
}
ARSimplex3D::
-ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices)
+ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices, const Delaunay& Dt)
{
Cell_handle c = f.first;
for (int i = 0; i < 4; ++i)
@@ -139,29 +135,33 @@
rho_ = squared_radius(p1, p2, p3);
attached_ = false;
- if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+ if (!Dt.is_infinite(c->vertex(f.second)) &&
+ CGAL::side_of_bounded_sphere(p1, p2, p3,
c->vertex(f.second)->point()) == CGAL::ON_BOUNDED_SIDE)
attached_ = true;
- else if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+ else if (!Dt.is_infinite(o->vertex(oi)) &&
+ CGAL::side_of_bounded_sphere(p1, p2, p3,
o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
attached_ = true;
else
alpha_ = rho_;
- SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
- SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
- if (c_iter == simplices.end()) // c is infinite
+ if (Dt.is_infinite(c))
{
+ SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
if (attached_) alpha_ = o_iter->first.alpha();
phi_const_ = o_iter->first.phi_const(); // FIXME: it's probably the other way around
}
- else if (o_iter == simplices.end()) // o is infinite
+ else if (Dt.is_infinite(o))
{
+ SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
if (attached_) alpha_ = c_iter->first.alpha();
phi_const_ = c_iter->first.phi_const(); // FIXME: it's probably the other way around
}
else
{
+ SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
+ SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
if (attached_) alpha_ = std::min(c_iter->first.alpha(), o_iter->first.alpha());
phi_const_ = std::min(c_iter->first.phi_const(), o_iter->first.phi_const());
}
@@ -241,10 +241,10 @@
rLog(rlARSimplex3D, "Vertices inserted");
for(Facet_iterator cur = Dt.finite_facets_begin(); cur != Dt.finite_facets_end(); ++cur)
- update_simplex_phi_map(ARSimplex3D(*cur, z, simplices), simplices);
+ update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt), simplices);
rLog(rlARSimplex3D, "Facets inserted");
for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
- update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt.incident_facets(*cur)), simplices);
+ update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt, Dt.incident_facets(*cur)), simplices);
rLog(rlARSimplex3D, "Edges inserted");
// Sort simplices by their alpha values
--- a/examples/ar-vineyard/ar-vineyard.cpp Wed Mar 19 12:50:35 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.cpp Wed Mar 19 13:50:51 2008 -0400
@@ -19,7 +19,7 @@
#endif
std::string infilename;
- double zx,zy,zz;
+ double zx,zy,zz,r;
std::string outfilename;
po::options_description hidden("Hidden options");
@@ -28,6 +28,7 @@
("x", po::value<double>(&zx), "x")
("y", po::value<double>(&zy), "y")
("z", po::value<double>(&zz), "z")
+ ("r", po::value<double>(&r), "r")
("output-file", po::value<std::string>(&outfilename), "Vineyard edges output");
std::vector<std::string> log_channels;
@@ -41,7 +42,7 @@
#endif
po::positional_options_description p;
- p.add("input-file", 1).add("x", 1).add("y", 1).add("z", 1).add("output-file", 1);
+ p.add("input-file", 1).add("x", 1).add("y", 1).add("z", 1).add("r", 1).add("output-file", 1);
po::options_description all; all.add(visible).add(hidden);
@@ -62,11 +63,12 @@
// Read command-line arguments
if (vm.count("help") || !vm.count("input-file") || !vm.count("output-file")
- || !vm.count("x") || !vm.count("y") || !vm.count("z"))
+ || !vm.count("x") || !vm.count("y") || !vm.count("z") || !vm.count("r"))
{
- std::cout << "Usage: ar-vineyard [OPTIONS] POINTS X Y Z OUTFILENAME" << std::endl;
+ std::cout << "Usage: ar-vineyard [OPTIONS] POINTS X Y Z R OUTFILENAME" << std::endl;
std::cout << " POINTS - filename containing points" << std::endl;
std::cout << " X,Y,Z - center-point z at which to compute the vineyard" << std::endl;
+ std::cout << " R - maximum radius" << std::endl;
std::cout << " OUTFILENAME - filename for the resulting vineyard" << std::endl;
std::cout << visible << std::endl;
std::cout << std::endl;
@@ -90,7 +92,7 @@
arv.compute_pairing();
// Compute vineyard
- arv.compute_vineyard();
+ arv.compute_vineyard(r);
std::cout << "Vineyard computed" << std::endl;
arv.vineyard()->save_edges(outfilename);
}
--- a/examples/ar-vineyard/ar-vineyard.h Wed Mar 19 12:50:35 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.h Wed Mar 19 13:50:51 2008 -0400
@@ -119,7 +119,7 @@
~ARVineyard();
void compute_pairing();
- void compute_vineyard();
+ void compute_vineyard(double max_radius);
const Filtration* filtration() const { return filtration_; }
const Vineyard* vineyard() const { return vineyard_; }
--- a/examples/ar-vineyard/ar-vineyard.hpp Wed Mar 19 12:50:35 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.hpp Wed Mar 19 13:50:51 2008 -0400
@@ -86,7 +86,7 @@
void
ARVineyard::
-compute_vineyard()
+compute_vineyard(double max_radius)
{
AssertMsg(filtration_->is_paired(), "Simplices must be paired for a vineyard to be computed");
@@ -113,7 +113,7 @@
// Simulate
change_evaluator(new KineticEvaluator(&simplex_sort_simulator, &trajectory_sort_simulator));
- while(!simplex_sort_simulator.reached_infinity() || !trajectory_sort_simulator.reached_infinity())
+ while ((simplex_sort_simulator.current_time() < max_radius || trajectory_sort_simulator.current_time() < max_radius) && !(simplex_sort_simulator.reached_infinity() && trajectory_sort_simulator.reached_infinity()))
{
if (*(simplex_sort_simulator.top()) < *(trajectory_sort_simulator.top()))
{