Fixed #c42 for ar-vineyard ar
authorDmitriy Morozov <morozov@cs.duke.edu>
Wed, 19 Mar 2008 13:50:51 -0400
branchar
changeset 86 73a54447b54a
parent 85 58deab1c8eb9
child 87 2c2e2f3b5d15
Fixed #c42 for ar-vineyard
.issues/c420501cc5285bbc
examples/ar-vineyard/ar-simplex3d.h
examples/ar-vineyard/ar-simplex3d.hpp
examples/ar-vineyard/ar-vineyard.cpp
examples/ar-vineyard/ar-vineyard.h
examples/ar-vineyard/ar-vineyard.hpp
--- 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()))
         {