diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_image_with_detection_of_features.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_detection_of_features.cpp index 7ef88f08410..c708bd25c43 100644 --- a/Mesh_3/examples/Mesh_3/mesh_3D_image_with_detection_of_features.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_3D_image_with_detection_of_features.cpp @@ -68,8 +68,7 @@ bool add_1D_features(const CGAL::Image_3& image, const double vx = image.vx(); const double vy = image.vy(); const double vz = image.vz(); - const double dist_bound = (std::min)(vx, - (std::min)(vy, vz)) / 256; + const double dist_bound = (std::min)(vx, (std::min)(vy, vz)) / 256; const double sq_dist_bound = dist_bound * dist_bound; const std::size_t xdim = image.xdim(); @@ -83,124 +82,124 @@ bool add_1D_features(const CGAL::Image_3& image, Del::Cell_handle start_cell; for (std::size_t k = 0, end_k = zdim - 1; k < end_k; ++k) - for (std::size_t j = 0, end_j = ydim - 1; j < end_j; ++j) - for (std::size_t i = 0, end_i = xdim - 1; i < end_i; ++i) - { - const K::Vector_3 translation{ i * vx, j * vy, k * vz }; - const Cube cube = { - static_evaluate(image.image(), i , j , k), - static_evaluate(image.image(), i + 1, j , k), - static_evaluate(image.image(), i , j + 1, k), - static_evaluate(image.image(), i + 1, j + 1, k), - static_evaluate(image.image(), i , j , k + 1), - static_evaluate(image.image(), i + 1, j , k + 1), - static_evaluate(image.image(), i , j + 1, k + 1), - static_evaluate(image.image(), i + 1, j + 1, k + 1), - }; /// TODO: optimize the access to the image data - bool monocolor = cube[0] == cube[1]; - for (int i = 2; i < 8; ++i) monocolor = monocolor && (cube[0] == cube[i]); - if (monocolor) continue; + for (std::size_t j = 0, end_j = ydim - 1; j < end_j; ++j) + for (std::size_t i = 0, end_i = xdim - 1; i < end_i; ++i) + { + const K::Vector_3 translation{ i * vx, j * vy, k * vz }; + const Cube cube = { + static_evaluate(image.image(), i , j , k), + static_evaluate(image.image(), i + 1, j , k), + static_evaluate(image.image(), i , j + 1, k), + static_evaluate(image.image(), i + 1, j + 1, k), + static_evaluate(image.image(), i , j , k + 1), + static_evaluate(image.image(), i + 1, j , k + 1), + static_evaluate(image.image(), i , j + 1, k + 1), + static_evaluate(image.image(), i + 1, j + 1, k + 1), + }; /// TODO: optimize the access to the image data + bool monocolor = cube[0] == cube[1]; + for (int i = 2; i < 8; ++i) monocolor = monocolor && (cube[0] == cube[i]); + if (monocolor) continue; - std::array inv_color_transformation{ INT_MIN, INT_MIN, INT_MIN, INT_MIN, - INT_MIN, INT_MIN, INT_MIN, INT_MIN }; - std::array color_transformation; - std::fill(color_transformation.begin(), color_transformation.end(), INT_MIN); - int nb_color = 0; - for (int i = 0; i < 8; ++i) { - if (color_transformation[cube[i]] == INT_MIN) { - color_transformation[cube[i]] = nb_color; - inv_color_transformation[nb_color] = cube[i]; - ++nb_color; - } - } - if (nb_color > 3) { - CGAL_warning_msg(nb_color > 3, "voxel with more than 3 colors"); - continue; - } - Cube reference_cube = { - (unsigned char)(color_transformation[cube[0]]), - (unsigned char)(color_transformation[cube[1]]), - (unsigned char)(color_transformation[cube[2]]), - (unsigned char)(color_transformation[cube[3]]), - (unsigned char)(color_transformation[cube[4]]), - (unsigned char)(color_transformation[cube[5]]), - (unsigned char)(color_transformation[cube[6]]), - (unsigned char)(color_transformation[cube[7]]), - }; - auto case_it = find_case(cases, reference_cube); - using std::end; - const bool case_found = (case_it != end(cases)); - if (case_found) reference_cube = combinations[(*case_it)[8]]; - else { -// std::cerr << "Warning: case not found: " << reference_cube << '\n'; - CGAL_error(); - }; -#ifdef CGAL_DEBUG_TRIPLE_LINES - std::cerr << "Cube " << cube << std::endl; - std::cerr << "reference cube " << reference_cube << std::endl; - std::cerr << " with transformation " << cube_isometries[(*case_it)[9]] << "\n"; -#endif // CGAL_DEBUG_TRIPLE_LINES - auto fct_it = lines.create_polylines_fcts.find(reference_cube); - if (fct_it != lines.create_polylines_fcts.end()) { -#ifdef CGAL_DEBUG_TRIPLE_LINES - std::cerr << "Using the function of " << Cube(fct_it->first) << "\n"; -#endif // CGAL_DEBUG_TRIPLE_LINES - Polylines cube_features = (fct_it->second)(10); - if (case_found) { - const Permutation& transformation = cube_isometries[(*case_it)[9]]; - Coordinates a1 = coordinates[transformation[0]]; - Coordinates u = coordinates[transformation[1]] - a1; - Coordinates v = coordinates[transformation[2]] - a1; - Coordinates w = coordinates[transformation[4]] - a1; - const Point_3 pa{ a1[0], a1[1], a1[2] }; - const Vector_3 vu{ u[0], u[1], u[2] }; - const Vector_3 vv{ v[0], v[1], v[2] }; - const Vector_3 vw{ w[0], w[1], w[2] }; -#ifdef CGAL_DEBUG_TRIPLE_LINES - std::cerr << "pa: " << pa << "\n"; - std::cerr << "vu: " << vu << "\n"; - std::cerr << "vv: " << vv << "\n"; - std::cerr << "vw: " << vw << "\n"; -#endif // CGAL_DEBUG_TRIPLE_LINES - for (auto& polyline : cube_features) { - for (auto& point : polyline) { - point = pa - + point.x() * vu - + point.y() * vv - + point.z() * vw; - point = { vx * point.x(), - vy * point.y(), - vz * point.z(), }; - point = point + translation; - } - for (int i = 0; i < 2; ++i) { - K::Point_3& extremity = (i == 0) ? polyline.front() : polyline.back(); - Del::Vertex_handle vh = triangulation.nearest_vertex(extremity, start_cell); - if (Del::Vertex_handle() != vh) { - if (squared_distance(vh->point(), extremity) < sq_dist_bound) { - extremity = vh->point(); - } - } - vh = triangulation.insert(extremity, start_cell); - start_cell = vh->cell(); - } - features_inside.push_back(std::move(polyline)); - } // end loop on polylines - } // end case where the transformation is not the identity - } // end if the reference_cube has polylines + std::array inv_color_transformation{ INT_MIN, INT_MIN, INT_MIN, INT_MIN, + INT_MIN, INT_MIN, INT_MIN, INT_MIN }; + std::array color_transformation; + std::fill(color_transformation.begin(), color_transformation.end(), INT_MIN); + int nb_color = 0; + for (int i = 0; i < 8; ++i) { + if (color_transformation[cube[i]] == INT_MIN) { + color_transformation[cube[i]] = nb_color; + inv_color_transformation[nb_color] = cube[i]; + ++nb_color; } + } + if (nb_color > 3) { + CGAL_warning_msg(nb_color > 3, "voxel with more than 3 colors"); + continue; + } + Cube reference_cube = { + (unsigned char)(color_transformation[cube[0]]), + (unsigned char)(color_transformation[cube[1]]), + (unsigned char)(color_transformation[cube[2]]), + (unsigned char)(color_transformation[cube[3]]), + (unsigned char)(color_transformation[cube[4]]), + (unsigned char)(color_transformation[cube[5]]), + (unsigned char)(color_transformation[cube[6]]), + (unsigned char)(color_transformation[cube[7]]), + }; + auto case_it = find_case(cases, reference_cube); + using std::end; + const bool case_found = (case_it != end(cases)); + if (case_found) reference_cube = combinations[(*case_it)[8]]; + else { + //std::cerr << "Warning: case not found: " << reference_cube << '\n'; + CGAL_error(); + }; +#ifdef CGAL_DEBUG_TRIPLE_LINES + std::cerr << "Cube " << cube << std::endl; + std::cerr << "reference cube " << reference_cube << std::endl; + std::cerr << " with transformation " << cube_isometries[(*case_it)[9]] << "\n"; +#endif // CGAL_DEBUG_TRIPLE_LINES + auto fct_it = lines.create_polylines_fcts.find(reference_cube); + if (fct_it != lines.create_polylines_fcts.end()) { +#ifdef CGAL_DEBUG_TRIPLE_LINES + std::cerr << "Using the function of " << Cube(fct_it->first) << "\n"; +#endif // CGAL_DEBUG_TRIPLE_LINES + Polylines cube_features = (fct_it->second)(10); + if (case_found) { + const Permutation& transformation = cube_isometries[(*case_it)[9]]; + Coordinates a1 = coordinates[transformation[0]]; + Coordinates u = coordinates[transformation[1]] - a1; + Coordinates v = coordinates[transformation[2]] - a1; + Coordinates w = coordinates[transformation[4]] - a1; + const Point_3 pa{ a1[0], a1[1], a1[2] }; + const Vector_3 vu{ u[0], u[1], u[2] }; + const Vector_3 vv{ v[0], v[1], v[2] }; + const Vector_3 vw{ w[0], w[1], w[2] }; +#ifdef CGAL_DEBUG_TRIPLE_LINES + std::cerr << "pa: " << pa << "\n"; + std::cerr << "vu: " << vu << "\n"; + std::cerr << "vv: " << vv << "\n"; + std::cerr << "vw: " << vw << "\n"; +#endif // CGAL_DEBUG_TRIPLE_LINES + for (auto& polyline : cube_features) { + for (auto& point : polyline) { + point = pa + + point.x() * vu + + point.y() * vv + + point.z() * vw; + point = { vx * point.x(), + vy * point.y(), + vz * point.z(), }; + point = point + translation; + } + for (int i = 0; i < 2; ++i) { + K::Point_3& extremity = (i == 0) ? polyline.front() : polyline.back(); + Del::Vertex_handle vh = triangulation.nearest_vertex(extremity, start_cell); + if (Del::Vertex_handle() != vh) { + if (squared_distance(vh->point(), extremity) < sq_dist_bound) { + extremity = vh->point(); + } + } + vh = triangulation.insert(extremity, start_cell); + start_cell = vh->cell(); + } + features_inside.push_back(std::move(polyline)); + } // end loop on polylines + } // end case where the transformation is not the identity + } // end if the reference_cube has polylines + } // call the split_graph_into_polylines, to create long polylines from the // short polylines that were generated per voxel. Polylines new_polylines_inside; CGAL::polylines_to_protect(new_polylines_inside, - features_inside.begin(), - features_inside.end()); + features_inside.begin(), + features_inside.end()); std::vector > polylines_on_bbox; CGAL::polylines_to_protect(image, polylines_on_bbox, - new_polylines_inside.begin(), - new_polylines_inside.end()); + new_polylines_inside.begin(), + new_polylines_inside.end()); domain.add_features(polylines_on_bbox.begin(), polylines_on_bbox.end()); @@ -211,9 +210,9 @@ bool add_1D_features(const CGAL::Image_3& image, std::ofstream output_polylines("out-generated.polylines.txt"); output_polylines.precision(17); for (auto poly : boost::range::join(polylines_on_bbox, new_polylines_inside)) { - output_polylines << poly.size(); - for (auto p : poly) output_polylines << " " << p; - output_polylines << std::endl; + output_polylines << poly.size(); + for (auto p : poly) output_polylines << " " << p; + output_polylines << std::endl; } return true;