refactored raycasting - inverted direction of raycasting - now each face is tested

fixed bug with custom seam drawings - square distance parameter named incorrectly

src/libslic3r/AABBTreeIndirect.hpp

@@ -759,8 +759,8 @@ inline bool is_any_triangle_in_radius(
         const TreeType 						&tree,
         // Point to which the closest point on the indexed triangle set is searched for.
         const VectorType					&point,
-        // Maximum distance in which triangle is search for
+        //Square of maximum distance in which triangle is searched for
-        typename VectorType::Scalar &max_distance)
+        typename VectorType::Scalar &max_distance_squared)
 {
     using Scalar = typename VectorType::Scalar;
     auto distancer = detail::IndexedTriangleSetDistancer<VertexType, IndexedFaceType, TreeType, VectorType>
@@ -774,7 +774,7 @@ inline bool is_any_triangle_in_radius(
 		return false;
 	}
 
-	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance, hit_idx, hit_point);
+	detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point);
 
     return hit_point.allFinite();
 }

src/libslic3r/GCode/SeamPlacerNG.cpp

@@ -162,6 +162,13 @@ Vec3f sample_sphere_uniform(const Vec2f &samples) {
         1.0f - 2.0f * samples.y()};
 }
 
+Vec3f sample_hemisphere_uniform(const Vec2f &samples) {
+    float term1 = 2.0f * float(PI) * samples.x();
+    float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y());
+    return {cos(term1) * term2, sin(term1) * term2,
+        abs(1.0f - 2.0f * samples.y())};
+}
+
 Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) {
     float term1 = 2.f * float(PI) * samples.x();
     float term2 = pow(samples.y(), 1.f / (power + 1.f));
@@ -170,88 +177,63 @@ Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) {
     return Vec3f(cos(term1) * term3, sin(term1) * term3, term2);
 }
 
-std::vector<HitInfo> raycast_visibility(const AABBTreeIndirect::Tree<3, float> &raycasting_tree,
+std::vector<FaceVisibilityInfo> raycast_visibility(const AABBTreeIndirect::Tree<3, float> &raycasting_tree,
         const indexed_triangle_set &triangles) {
-
-    auto bbox = raycasting_tree.node(0).bbox;
-    Vec3f vision_sphere_center = bbox.center().cast<float>();
-    Vec3f side_sizes = bbox.sizes().cast<float>();
-    float vision_sphere_raidus = (side_sizes.norm() * 0.55); // 0.5 (half) covers whole object,
-    // 0.05 added to avoid corner cases
-    // very rough approximation of object surface area from its bounding sphere
-    float approx_area = 4 * PI * vision_sphere_raidus * vision_sphere_raidus;
-    float local_considered_area = PI * SeamPlacer::considered_area_radius * SeamPlacer::considered_area_radius;
-    size_t ray_count = SeamPlacer::expected_hits_per_area * approx_area / local_considered_area;
-
-    // Prepare random samples per ray
-    // std::random_device rnd_device;
-    // use fixed seed, we can backtrack potential issues easier
-    std::mt19937 mersenne_engine { 131 };
-    std::uniform_real_distribution<float> dist { 0, 1 };
-
-    auto gen = [&dist, &mersenne_engine]() {
-        return Vec2f(dist(mersenne_engine), dist(mersenne_engine));
-    };
-
     BOOST_LOG_TRIVIAL(debug)
-    << "SeamPlacer: generate random samples: start";
+    << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: start";
-    std::vector<Vec2f> global_dir_random_samples(ray_count);
-    generate(begin(global_dir_random_samples), end(global_dir_random_samples), gen);
-    std::vector<Vec2f> local_dir_random_samples(ray_count);
-    generate(begin(local_dir_random_samples), end(local_dir_random_samples), gen);
-    BOOST_LOG_TRIVIAL(debug)
-    << "SeamPlacer: generate random samples: end";
 
-    BOOST_LOG_TRIVIAL(debug)
+    float step_size = 1.0f / SeamPlacer::sqr_rays_per_triangle;
-    << "SeamPlacer: raycast visibility for " << ray_count << " rays: start";
+    std::vector<Vec3f> precomputed_sample_directions(
+            SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle);
+    for (size_t x_idx = 0; x_idx < SeamPlacer::sqr_rays_per_triangle; ++x_idx) {
+        float sample_x = x_idx * step_size + step_size / 2.0;
+        for (size_t y_idx = 0; y_idx < SeamPlacer::sqr_rays_per_triangle; ++y_idx) {
+            size_t dir_index = x_idx * SeamPlacer::sqr_rays_per_triangle + y_idx;
+            float sample_y = y_idx * step_size + step_size / 2.0;
+            precomputed_sample_directions[dir_index] = sample_hemisphere_uniform( { sample_x, sample_y });
+        }
+    }
 
-    std::vector<HitInfo> hit_points = tbb::parallel_reduce(tbb::blocked_range<size_t>(0, ray_count),
+    std::vector<FaceVisibilityInfo> result(triangles.indices.size());
-            std::vector<HitInfo> { },
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, result.size()),
-            [&](tbb::blocked_range<size_t> r, std::vector<HitInfo> init) {
+            [&](tbb::blocked_range<size_t> r) {
-                for (size_t index = r.begin(); index < r.end(); ++index) {
+                for (size_t face_index = r.begin(); face_index < r.end(); ++face_index) {
-                    //generate global ray direction
+                    FaceVisibilityInfo &dest = result[face_index];
-                    Vec3f global_ray_dir = sample_sphere_uniform(global_dir_random_samples[index]);
+                    dest.visibility = 1.0f;
-                    //place the ray origin on the bounding sphere
+                    constexpr float decrease =  1.0f / (SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle);
-                    Vec3f ray_origin = (vision_sphere_center - global_ray_dir * vision_sphere_raidus);
-                    // compute local ray direction as cosine hemisphere sample - the rays dont aim directly to the middle
-                    Vec3f local_dir = sample_power_cosine_hemisphere(local_dir_random_samples[index], SeamPlacer::cosine_hemisphere_sampling_power);
 
+
+                    Vec3i face = triangles.indices[face_index];
+                    Vec3f A = triangles.vertices[face.x()];
+                    Vec3f B = triangles.vertices[face.y()];
+                    Vec3f C = triangles.vertices[face.z()];
+                    Vec3f center = (A + B + C) / 3.0f;
+                    Vec3f normal = ((B - A).cross(C - B)).normalized();
                     // apply the local direction via Frame struct - the local_dir is with respect to +Z being forward
                     Frame f;
-                    f.set_from_z(global_ray_dir);
+                    f.set_from_z(normal);
-                    Vec3f final_ray_dir = (f.to_world(local_dir));
 
+                    for (const auto &dir : precomputed_sample_directions) {
+                        Vec3f final_ray_dir = (f.to_world(dir));
                         igl::Hit hitpoint;
                         // FIXME: This AABBTTreeIndirect query will not compile for float ray origin and
                         // direction.
-                    Vec3d ray_origin_d = ray_origin.cast<double>();
+                        Vec3d ray_origin_d = (center+0.1*normal).cast<double>();
                         Vec3d final_ray_dir_d = final_ray_dir.cast<double>();
                         bool hit = AABBTreeIndirect::intersect_ray_first_hit(triangles.vertices,
                                 triangles.indices, raycasting_tree, ray_origin_d, final_ray_dir_d, hitpoint);
 
                         if (hit) {
-                        auto face = triangles.indices[hitpoint.id];
+                            dest.visibility -= decrease;
-                        auto edge1 = triangles.vertices[face[1]] - triangles.vertices[face[0]];
-                        auto edge2 = triangles.vertices[face[2]] - triangles.vertices[face[0]];
-
-                        Vec3f hit_pos = (triangles.vertices[face[0]] + edge1 * hitpoint.u + edge2 * hitpoint.v);
-                        Vec3f surface_normal = its_face_normal(triangles, hitpoint.id);
-
-                        init.push_back(HitInfo { hit_pos, surface_normal });
                         }
                     }
-                return init;
-            },
-            [](std::vector<HitInfo> left, const std::vector<HitInfo>& right) {
-                left.insert(left.end(), right.begin(), right.end());
-                return left;
                 }
-            );
+            });
 
     BOOST_LOG_TRIVIAL(debug)
-    << "SeamPlacer: raycast visibility for " << ray_count << " rays: end";
+    << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: end";
 
-    return hit_points;
+    return result;
 }
 
 std::vector<float> calculate_polygon_angles_at_vertices(const Polygon &polygon, const std::vector<float> &lengths,
@@ -318,54 +300,41 @@ std::vector<float> calculate_polygon_angles_at_vertices(const Polygon &polygon,
 
 // structure to store global information about the model - occlusion hits, enforcers, blockers
 struct GlobalModelInfo {
-    std::vector<HitInfo> geometry_raycast_hits;
+    indexed_triangle_set model;
-    KDTreeIndirect<3, float, HitInfoCoordinateFunctor> raycast_hits_tree;
+    AABBTreeIndirect::Tree<3, float> model_tree;
+    std::vector<FaceVisibilityInfo> visiblity_info;
     indexed_triangle_set enforcers;
     indexed_triangle_set blockers;
     AABBTreeIndirect::Tree<3, float> enforcers_tree;
     AABBTreeIndirect::Tree<3, float> blockers_tree;
 
-    GlobalModelInfo() :
-            raycast_hits_tree(HitInfoCoordinateFunctor { &geometry_raycast_hits }) {
-    }
-
     bool is_enforced(const Vec3f &position, float radius) const {
         if (enforcers.empty()) {
             return false;
         }
+        float radius_sqr = radius*radius;
         return AABBTreeIndirect::is_any_triangle_in_radius(enforcers.vertices, enforcers.indices,
-                enforcers_tree, position, radius);
+                enforcers_tree, position, radius_sqr);
     }
 
     bool is_blocked(const Vec3f &position, float radius) const {
         if (blockers.empty()) {
             return false;
         }
+        float radius_sqr = radius*radius;
         return AABBTreeIndirect::is_any_triangle_in_radius(blockers.vertices, blockers.indices,
-                blockers_tree, position, radius);
+                blockers_tree, position, radius_sqr);
     }
 
     float calculate_point_visibility(const Vec3f &position) const {
-        size_t closest_point_index = find_closest_point(raycast_hits_tree, position);
+        size_t hit_idx;
-        if (closest_point_index == raycast_hits_tree.npos
+        Vec3f hit_point;
-                ||
+        if (AABBTreeIndirect::squared_distance_to_indexed_triangle_set(model.vertices, model.indices, model_tree,
-                (position - geometry_raycast_hits[closest_point_index].position).norm()
+                position, hit_idx, hit_point) >= 0) {
-                        > SeamPlacer::considered_area_radius) {
+            return visiblity_info[hit_idx].visibility;
-            return 0;
+        } else {
+            return 0.0f;
         }
-        auto nearby_points = find_nearby_points(raycast_hits_tree, position, SeamPlacer::considered_area_radius);
-        Vec3f local_normal = geometry_raycast_hits[closest_point_index].surface_normal;
-
-        float visibility = 0;
-        for (const auto &hit_point_index : nearby_points) {
-            float dot_product = local_normal.dot(geometry_raycast_hits[hit_point_index].surface_normal);
-            if (dot_product > 0) {
-                float distance =
-                        (position - geometry_raycast_hits[hit_point_index].position).norm();
-                visibility += dot_product*gauss(distance, 0.0f, 1.0f, 0.5f);
-            }
-        }
-        return visibility;
 
     }
 
@@ -373,7 +342,7 @@ struct GlobalModelInfo {
     void debug_export(const indexed_triangle_set &obj_mesh, const char *file_name) const {
         indexed_triangle_set divided_mesh = obj_mesh;       //  subdivide(obj_mesh, SeamPlacer::considered_area_radius);
         Slic3r::CNumericLocalesSetter locales_setter;
-        {
+
         FILE *fp = boost::nowide::fopen(file_name, "w");
         if (fp == nullptr) {
             BOOST_LOG_TRIVIAL(error)
@@ -383,7 +352,8 @@ struct GlobalModelInfo {
 
         for (size_t i = 0; i < divided_mesh.vertices.size(); ++i) {
             float visibility = calculate_point_visibility(divided_mesh.vertices[i]);
-                Vec3f color =  value_to_rgbf(0, SeamPlacer::expected_hits_per_area, visibility);
+            Vec3f color = value_to_rgbf(0.0f, 1.0f,
+                    visibility);
             fprintf(fp, "v %f %f %f  %f %f %f\n",
                     divided_mesh.vertices[i](0), divided_mesh.vertices[i](1), divided_mesh.vertices[i](2),
                     color(0), color(1), color(2)
@@ -393,24 +363,7 @@ struct GlobalModelInfo {
             fprintf(fp, "f %d %d %d\n", divided_mesh.indices[i][0] + 1, divided_mesh.indices[i][1] + 1,
                     divided_mesh.indices[i][2] + 1);
         fclose(fp);
-        }
-        {
-            auto fname = std::string("hits_").append(file_name);
-            FILE *fp = boost::nowide::fopen(fname.c_str(), "w");
-            if (fp == nullptr) {
-                BOOST_LOG_TRIVIAL(error)
-                << "Couldn't open " << fname << " for writing";
-            }
 
-            for (size_t i = 0; i < geometry_raycast_hits.size(); ++i) {
-                Vec3f surface_normal = (geometry_raycast_hits[i].surface_normal + Vec3f(1.0, 1.0, 1.0)) / 2.0;
-                fprintf(fp, "v %f %f %f %f %f %f \n", geometry_raycast_hits[i].position[0],
-                        geometry_raycast_hits[i].position[1],
-                        geometry_raycast_hits[i].position[2], surface_normal[0], surface_normal[1],
-                        surface_normal[2]);
-            }
-            fclose(fp);
-        }
     }
 #endif
 }
@@ -512,7 +465,7 @@ void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, std::
             if (global_model_info.is_enforced(position, distance_to_next)
                     || global_model_info.is_blocked(position, distance_to_next)) {
                 Vec3f vec_to_next = (pos_of_next - position).normalized();
-                float step_size = SeamPlacer::enforcer_blocker_distance_tolerance / 2.0f;
+                float step_size = SeamPlacer::enforcer_blocker_oversampling_distance;
                 float step = step_size;
                 while (step < distance_to_next) {
                     oversampled_points.push(position + vec_to_next * step);
@@ -620,16 +573,17 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) {
             its_merge(triangle_set, model_its);
         }
     }
-
-    its_transform(triangle_set, obj_transform);
     float target_error = SeamPlacer::raycasting_decimation_target_error;
     its_quadric_edge_collapse(triangle_set, 0, &target_error, nullptr, nullptr);
+    triangle_set = subdivide(triangle_set, SeamPlacer::raycasting_subdivision_target_length);
+    its_transform(triangle_set, obj_transform);
 
     auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices,
             triangle_set.indices);
 
-    result.geometry_raycast_hits = raycast_visibility(raycasting_tree, triangle_set);
+    result.model = triangle_set;
-    result.raycast_hits_tree.build(result.geometry_raycast_hits.size());
+    result.model_tree = raycasting_tree;
+    result.visiblity_info = raycast_visibility(raycasting_tree, triangle_set);
 
     BOOST_LOG_TRIVIAL(debug)
     << "SeamPlacer: build AABB tree for raycasting and gather occlusion info: end";
@@ -679,9 +633,8 @@ struct SeamComparator {
             setup(setup) {
     }
 
-
     float compute_angle_penalty(float ccw_angle) const {
-        return gauss(ccw_angle, 0.2f, 1.0f, 0.7f);
+        return gauss(ccw_angle, 0.2f, 1.0f, 4.0f);
     }
 
     // Standard comparator, must respect the requirements of comparators (e.g. give same result on same inputs) for sorting usage
@@ -704,8 +657,8 @@ struct SeamComparator {
             return a.position.y() > b.position.y();
         }
 
-        return (a.visibility + SeamPlacer::expected_hits_per_area) * compute_angle_penalty(a.local_ccw_angle) <
+        return (a.visibility + 1.0f) * compute_angle_penalty(a.local_ccw_angle) <
-                (b.visibility + SeamPlacer::expected_hits_per_area) * compute_angle_penalty(b.local_ccw_angle);
+                (b.visibility + 1.0f) * compute_angle_penalty(b.local_ccw_angle);
     }
 
     // Comparator used during alignment. If there is close potential aligned point, it is comapred to the current
@@ -728,17 +681,17 @@ struct SeamComparator {
             return a.position.y() > b.position.y();
         }
 
-        return (a.visibility + SeamPlacer::expected_hits_per_area) * compute_angle_penalty(a.local_ccw_angle) * 0.8f <=
+        return (a.visibility +  1.0f) * compute_angle_penalty(a.local_ccw_angle) * 0.5f <=
-                (b.visibility + SeamPlacer::expected_hits_per_area) * compute_angle_penalty(b.local_ccw_angle);
+                (b.visibility +  1.0f) * compute_angle_penalty(b.local_ccw_angle);
     }
 
+    //returns negative value of penalties, should be nromalized against others in the same perimeter for use
     float get_weight(const SeamCandidate &a) const {
         if (setup == SeamPosition::spRear) {
             return a.position.y();
         }
 
-        //return negative, beacuse we want to minimize the absolute of this value (sort of antiweight). normalization in alignment fixes that with respect to other points
+        return -(a.visibility + 1.0f) * compute_angle_penalty(a.local_ccw_angle);
-        return -(a.visibility + SeamPlacer::expected_hits_per_area) * compute_angle_penalty(a.local_ccw_angle);
     }
 }
 ;
@@ -776,7 +729,8 @@ void debug_export_points(const std::vector<std::vector<SeamPlacerImpl::SeamCandi
             visibility_svg.draw(scaled(Vec2f(point.position.head<2>())), visibility_fill);
 
             Vec3i weight_color = value_rgbi(min_weight, max_weight, comparator.get_weight(point));
-            std::string weight_fill = "rgb(" + std::to_string(weight_color.x()) + "," + std::to_string(weight_color.y()) + ","
+            std::string weight_fill = "rgb(" + std::to_string(weight_color.x()) + "," + std::to_string(weight_color.y())
+                    + ","
                     + std::to_string(weight_color.z()) + ")";
             weight_svg.draw(scaled(Vec2f(point.position.head<2>())), weight_fill);
         }
@@ -1181,7 +1135,8 @@ void SeamPlacer::init(const Print &print) {
             << "SeamPlacer: align_seam_points : end";
         }
 
-        debug_export_points(m_perimeter_points_per_object[po], po->bounding_box(), std::to_string(po->id().id), comparator);
+        debug_export_points(m_perimeter_points_per_object[po], po->bounding_box(), std::to_string(po->id().id),
+                comparator);
     }
 }
 

src/libslic3r/GCode/SeamPlacerNG.hpp

@@ -65,10 +65,8 @@ struct SeamCandidate {
     EnforcedBlockedSeamPoint type;
 };
 
-// struct to represent hits of the mesh during occulision raycasting.
+struct FaceVisibilityInfo {
-struct HitInfo {
+    float visibility;
-    Vec3f position;
-    Vec3f surface_normal;
 };
 
 struct SeamCandidateCoordinateFunctor {
@@ -80,38 +78,24 @@ struct SeamCandidateCoordinateFunctor {
         return seam_candidates->operator[](index).position[dim];
     }
 };
-
-struct HitInfoCoordinateFunctor {
-    HitInfoCoordinateFunctor(std::vector<HitInfo> *hit_points) :
-            hit_points(hit_points) {
-    }
-    std::vector<HitInfo> *hit_points;
-    float operator()(size_t index, size_t dim) const {
-        return hit_points->operator[](index).position[dim];
-    }
-};
 } // namespace SeamPlacerImpl
 
 class SeamPlacer {
 public:
     using SeamCandidatesTree =
     KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>;
-    // Rough estimates of hits of the mesh during raycasting per surface circle defined by considered_area_radius
+    static constexpr float raycasting_decimation_target_error = 2.0f;
-    static constexpr float expected_hits_per_area = 1000.0f;
+    static constexpr float raycasting_subdivision_target_length = 3.0f;
-    // area considered when computing number of rays and then gathering visiblity info from the hits
+    //square of number of rays per triangle
-    static constexpr float considered_area_radius = 3.0f;
+    static constexpr size_t sqr_rays_per_triangle = 8;
-    // quadric error limit of quadric decimation function used on the mesh before raycasting
-    static constexpr float raycasting_decimation_target_error = 0.5f;
-
-    // cosine sampling power represents how prefered are forward directions when raycasting from given spot
-    // in this case, forward direction means towards the center of the mesh
-    static constexpr float cosine_hemisphere_sampling_power = 5.0f;
 
     // arm length used during angles computation
     static constexpr float polygon_local_angles_arm_distance = 0.4f;
 
     // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer
-    static constexpr float enforcer_blocker_distance_tolerance = 0.2f;
+    static constexpr float enforcer_blocker_distance_tolerance = 0.3f;
+    // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size
+    static constexpr float enforcer_blocker_oversampling_distance = 0.3f;
 
     // When searching for seam clusters for alignment:
     // seam_align_tolerable_dist - if seam is closer to the previous seam position projected to the current layer than this value,