From 1d5e9e915079ef33d309d8ad0b90187bbb18fb6f Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 15 Jan 2018 16:14:26 +0100 Subject: [PATCH] WIP: implementation for scalar image ... tested only on segmented images! --- .../CGAL/Mesh_3/polylines_to_protect.h | 312 +++++++++++++----- 1 file changed, 229 insertions(+), 83 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h index 1fc27924116..26963b42b80 100644 --- a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h +++ b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h @@ -37,6 +37,7 @@ // CGAL::Null_subdomain_index #include // for boost::prior #include +#include #include #include @@ -58,6 +59,116 @@ struct Returns_midpoint { = K().construct_midpoint_3_object(); return a < b ? midpt(a, b) : midpt(b, a); } + double operator()(const NT, + const NT) const + { + return 0.5; + } +}; + +class Isoline_equation { + double a; + double b; + double c; + double d; + +public: + Isoline_equation(double a, double b, double c, double d) + : a(a) + , b(b) + , c(c) + , d(d) + {} + + double y(double x) const { + return ( x*(a-b)-a ) / + ( c-a+x*(a-b-c+d) ); + // If (a+d) == (b+c), then the curve is a straight line. + } + double x(double y) const { + return ( y*(a-c)-a ) / + ( b-a+y*(a-b-c+d) ); + // If (a+d) == (b+c), then the curve is a straight line. + } +}; // end of class Isoline_equation + +template +class Insert_curve_in_graph { + typedef typename Geom_traits::Point_3 Point; + typedef typename Geom_traits::Vector_3 Vector; + + G_manip& g_manip; + const Geom_traits& gt; + const typename Geom_traits::Construct_translated_point_3 translate; + const typename Geom_traits::Construct_scaled_vector_3 scale; + const int prec; + +public: + Insert_curve_in_graph(G_manip& g_manip, + const Geom_traits& gt, + const int prec = 10) + : g_manip(g_manip) + , gt(gt) + , translate(gt.construct_translated_point_3_object()) + , scale(gt.construct_scaled_vector_3_object()) + , prec(prec) + {} + + struct Coords { + double x; + double y; + }; + + void insert_curve(Isoline_equation equation, + vertex_descriptor start_v, + vertex_descriptor end_v, + Coords start, + Coords end, + Point p00, + Vector vx, + Vector vy) + { + if(CGAL::abs(start.x - end.x) >= CGAL::abs(start.y - end.y)) { + insert_curve(equation, &Isoline_equation::y, + start_v, end_v, + start.x, end.x, p00, vx, vy); + } else { + insert_curve(equation, &Isoline_equation::x, + start_v, end_v, + start.y, end.y, p00, vy, vx); + } + } +private: + void insert_curve(Isoline_equation equation, + double (Isoline_equation::* f)(double) const, + vertex_descriptor start_v, + vertex_descriptor end_v, + double start, + double end, + Point p00, + Vector vx, + Vector vy) + { + const double step = (end - start)/prec; + const double stop = end-step/2; + const bool step_is_positive = (step > 0); + vertex_descriptor old = start_v; + vertex_descriptor v_int; + for(double x = start + step; + (step_is_positive ? x < stop : x > stop); + x += step) + { + const double y = (equation.*f)(x); + const Point inter_p = + translate(translate(p00, + scale(vx, x)), + scale(vy, y)); + v_int = g_manip.get_vertex(inter_p); + g_manip.try_add_edge(old, v_int); + old = v_int; + } + g_manip.try_add_edge(v_int, end_v); + } }; template void -polylines_to_protect(const CGAL::Image_3& cgal_image, - const double vx, const double vy, const double vz, - std::vector >& polylines, - Image_word_type*, - Null_subdomain_index null, - DomainFunctor domain_fct, - InterpolationFunctor interpolate, - PolylineInputIterator existing_polylines_begin, - PolylineInputIterator existing_polylines_end) +polylines_to_protect +(const CGAL::Image_3& cgal_image, + const double vx, const double vy, const double vz, + std::vector >& polylines, + Image_word_type*, + Null_subdomain_index null, + DomainFunctor domain_fct, + InterpolationFunctor interpolate, + PolylineInputIterator existing_polylines_begin, + PolylineInputIterator existing_polylines_end, + boost::optional scalar_interpolation_value = boost::none, + int prec = 10) { typedef typename DomainFunctor::result_type Domain_type; typedef typename Kernel_traits

::Kernel K; @@ -249,10 +363,18 @@ polylines_to_protect(const CGAL::Image_3& cgal_image, const int image_dims[3] = { xdim, ydim, zdim }; Graph graph; - internal::Mesh_3::Graph_manipulations g_manip(graph, interpolate); + + using namespace CGAL::internal::Mesh_3; + + typedef internal::Mesh_3::Graph_manipulations G_manip; + G_manip g_manip(graph, interpolate); + + Insert_curve_in_graph insert_curve_in_graph(g_manip, K(), prec); std::size_t case4 = 0, // 4 colors @@ -306,9 +428,8 @@ polylines_to_protect(const CGAL::Image_3& cgal_image, { pix11, Point_3(), Domain_type(), 0 } }} }}; std::map pixel_values_set; - for(int ii = 0; ii < 2; ++ii) - for(int jj = 0; jj < 2; ++jj) - { + for(int ii = 0; ii < 2; ++ii) { + for(int jj = 0; jj < 2; ++jj) { const Pixel& pixel = square[ii][jj].pixel; double x = pixel[0] * vx; double y = pixel[1] * vy; @@ -319,8 +440,15 @@ polylines_to_protect(const CGAL::Image_3& cgal_image, pixel[1], pixel[2])); square[ii][jj].domain = domain_fct(square[ii][jj].word); + if(scalar_interpolation_value != boost::none) { + square[ii][jj].word = + Image_word_type(square[ii][jj].word - + (*scalar_interpolation_value)); + } ++pixel_values_set[square[ii][jj].domain]; } + } + const Point_3& p00 = square[0][0].point; const Point_3& p10 = square[1][0].point; const Point_3& p01 = square[0][1].point; @@ -375,6 +503,7 @@ polylines_to_protect(const CGAL::Image_3& cgal_image, CGAL_assertion(square[1][0].domain != square[1][1].domain); CGAL_assertion(square[0][1].domain != square[1][1].domain); CGAL_assertion(square[0][1].domain != square[1][0].domain); + case_4: // case 4 or case 2-2 // @@ -393,7 +522,9 @@ case_4: vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11); vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10); - vertex_descriptor vmid = g_manip.get_vertex(midpoint(p00, p11)); + vertex_descriptor vmid = g_manip.split(graph[left], graph[right], + v00, v10, + false, false); g_manip.try_add_edge(left , vmid); g_manip.try_add_edge(right , vmid); g_manip.try_add_edge(top , vmid); @@ -432,37 +563,34 @@ case_4: CGAL_assertion(square[1][1].domain != square[0][0].domain); CGAL_assertion(square[0][1].domain != square[0][0].domain); CGAL_assertion(square[0][1].domain != square[1][1].domain); +case_1_2_1: ++case121; vertex_descriptor left = g_manip.split(p00, p01, v00, v01, out00, out01); vertex_descriptor right = g_manip.split(p10, p11, v10, v11, out10, out11); vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11); vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10); - vertex_descriptor old_left = left; - vertex_descriptor old_right = right; - vertex_descriptor v_int_left, v_int_right; - - // approximate the arcs by 10 segments - // -> 9 intermediate vertices - for(double x = 0.05; x < 0.5; x+= 0.05) - { - const Point_3 inter_left = - translate(p00 - , x * vector(p00, p10) // x - + ((1.-2.*x)/(2.-3.*x)) * vector(p00, p01)); // y - const Point_3 inter_right = - translate(p11 - , x * vector(p11, p01) // x - + ((1.-2.*x)/(2.-3.*x)) * vector(p11, p10)); // y - v_int_left = g_manip.get_vertex(inter_left); - v_int_right = g_manip.get_vertex(inter_right); - g_manip.try_add_edge(old_left, v_int_left); - g_manip.try_add_edge(old_right, v_int_right); - old_left = v_int_left; - old_right = v_int_right; + Isoline_equation equation = + (scalar_interpolation_value == boost::none) ? + Isoline_equation(1, -1, -1, 0) : + Isoline_equation(v00, v10, v01, v11); + insert_curve_in_graph.insert_curve(equation, + left, bottom, + {0., interpolate(v00, v01) }, + {interpolate(v00, v10), 0. }, + p00, + p10 - p00, + p01 - p00); + if(scalar_interpolation_value == boost::none) { + equation = Isoline_equation(0, -1, -1, 1); } - g_manip.try_add_edge(v_int_left, bottom); - g_manip.try_add_edge(v_int_right, top); + insert_curve_in_graph.insert_curve(equation, + top, right, + {interpolate(v01, v11), 1. }, + {1., interpolate(v10, v11) }, + p00, + p10 - p00, + p01 - p00); } else { // case 2-1-1 if(square[0][0].domain == square[1][0].domain) { @@ -483,6 +611,8 @@ case_4: CGAL_assertion(square[0][0].domain != square[1][1].domain); CGAL_assertion(square[1][0].domain != square[1][1].domain); ++case211; + // Note: this case cannot occur with non-segmented scalar + // images, because it needs three domains. // // A-------------B // | | | @@ -514,32 +644,20 @@ case_4: vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11); vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10); - vertex_descriptor old_top = top; - vertex_descriptor old_bottom = bottom; - vertex_descriptor v_int_top, v_int_bottom; - - // approximate the arcs by 10 segments - // -> 9 intermediate vertices - for(double x = 0.51666; x < 0.66; x+= 0.016666) - { - const Point_3 inter_top = - translate(p00 - , x * vector(p00, p10) // x - + ((1./x) - 1.) * vector(p00, p01)); // y - const Point_3 inter_bottom = - translate(p00 - , x * vector(p00, p10) // x - + (2.-(1./x)) * vector(p00, p01)); // y - v_int_top = g_manip.get_vertex(inter_top); - v_int_bottom = g_manip.get_vertex(inter_bottom); - g_manip.try_add_edge(old_top, v_int_top); - g_manip.try_add_edge(old_bottom, v_int_bottom); - old_top = v_int_top; - old_bottom = v_int_bottom; - } - - g_manip.try_add_edge(v_int_bottom, v_inter); - g_manip.try_add_edge(v_int_top, v_inter); + insert_curve_in_graph.insert_curve(Isoline_equation(1, -1, 1, 0), + bottom, v_inter, + { 0.5 , 0. }, + { 2./3., 0.5 }, + p00, + p10 - p00, + p01 - p00); + insert_curve_in_graph.insert_curve(Isoline_equation(1, 0, 1, -1), + top, v_inter, + { 0.5 , 1. }, + { 2./3., 0.5 }, + p00, + p10 - p00, + p01 - p00); g_manip.try_add_edge(right, v_inter); } // end case 2-1-1 } // end `case 3:` @@ -585,13 +703,43 @@ case_4: CGAL_assertion(square[1][0].domain!=square[0][1].domain); vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11); vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10); - g_manip.try_add_edge(top, bottom); + if(scalar_interpolation_value == boost::none) { + g_manip.try_add_edge(top, bottom); + } else { + insert_curve_in_graph.insert_curve(Isoline_equation(v00, v10, + v01, v11), + top, bottom, + {interpolate(v01,v11), 1.}, + {interpolate(v00,v10), 0.}, + p00, + p10 - p00, + p01 - p00); + } } else { // Else diagonal case case 2-2 // Same as the case with 4 colors CGAL_assertion(square[0][0].domain==square[1][1].domain); CGAL_assertion(square[1][0].domain==square[0][1].domain); CGAL_assertion(square[0][0].domain!=square[0][1].domain); + + if(scalar_interpolation_value != boost::none) { + // Compute the squared distance between the two branches of + // the hyperbola. + const double discrimant = double(v00) * v11 - double(v01) * v10; + const double squared_distance = + 8 * CGAL::abs(discrimant) / CGAL::square(double(v00) - v10 - v01 + v11); + if(CGAL::square(prec) * squared_distance > 1) + { + // In that case, the case 1-2-1 will be applied + if(discrimant > 0.) { + // Vertical swap + std::swap(square[0][1], square[0][0]); std::swap(out01, out00); + std::swap(square[1][1], square[1][0]); std::swap(out11, out10); + } + goto case_1_2_1; + } + } // scalar interpolation + goto case_4; } } @@ -639,22 +787,20 @@ case_4: CGAL_assertion(square[1][1].domain == square[0][0].domain); CGAL_assertion(square[0][1].domain == square[0][0].domain); vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10); - vertex_descriptor old = bottom; - - vertex_descriptor v_int; - for(double x = 0.55; x < 1.; x+= 0.05) - { - const Point_3 inter = - translate(p00 - , x * vector(p00, p10) // x - + (1.-1./(2.*x)) * vector(p00, p01)); // y - v_int = g_manip.get_vertex(inter); - g_manip.try_add_edge(old, v_int); - old = v_int; - } - vertex_descriptor right = g_manip.split(p10, p11, v10, v11, out10, out11); - g_manip.try_add_edge(v_int, right); + + Isoline_equation equation = + (scalar_interpolation_value == boost::none) ? + Isoline_equation(1, -1, 1, 1) : + Isoline_equation(v00, v10, v01, v11); + + insert_curve_in_graph.insert_curve(equation, + bottom, right, + {interpolate(v00, v10), 0. }, + {1., interpolate(v10, v11) }, + p00, + p10 - p00, + p01 - p00); } } break;