C++ MatrixFr::rows方法代码示例

void SimpleInflator::generate_joint(
        const MatrixFr& pts, const VectorI& source_ids, size_t vertex_index) {
    const size_t dim = m_wire_network->get_dim();
    ConvexHullEngine::Ptr convex_hull = ConvexHullEngine::create(dim, "qhull");
    convex_hull->run(pts);

    MatrixFr vertices = convex_hull->get_vertices();
    MatrixIr faces = convex_hull->get_faces();
    VectorI  index_map = convex_hull->get_index_map();

    if (dim == 2) {
        // Need to triangulate the loop.
        const size_t num_vertices = vertices.rows();
        TriangleWrapper tri(vertices, faces);
        tri.run(std::numeric_limits<Float>::max(), false, false, false);
        vertices = tri.get_vertices();
        faces = tri.get_faces();
        assert(vertices.rows() == num_vertices);
    }

    m_vertex_list.push_back(vertices);
    const size_t num_faces = faces.rows();
    for (size_t i=0; i<num_faces; i++) {
        const auto& face = faces.row(i);
        if (dim == 3) {
            auto ori_indices = map_indices(face, index_map);
            if (belong_to_the_same_loop(ori_indices, source_ids)) continue;
        }
        m_face_list.push_back(face.array() + m_num_vertex_accumulated);
        m_face_source_list.push_back(vertex_index+1);
    }

    m_num_vertex_accumulated += vertices.rows();
}

void GeometryCorrectionTable::apply_z_correction(
        const Vector3F& edge_dir, MatrixFr& loop) {
    //const Float max_z_error = 0.125;
    //const Float max_z_error = 0.09;
    const Float max_z_error = 0.00;
    VectorF bbox_min = loop.colwise().minCoeff();
    VectorF bbox_max = loop.colwise().maxCoeff();
    VectorF bbox_center = 0.5 * (bbox_min + bbox_max);

    Vector3F side_dir = edge_dir.cross(Vector3F::UnitZ());
    Float sin_val = side_dir.norm();
    if (sin_val < 1e-3) return;

    const size_t num_vts = loop.rows();
    for (size_t i=0; i<num_vts; i++) {
        Vector3F v = loop.row(i) - bbox_center.transpose();
        Float side_component = side_dir.dot(v) / sin_val;
        Vector3F proj_v = v - side_component * side_dir / sin_val;
        Float proj_component = proj_v.norm();
        if (proj_component > 1e-3) {
            proj_v -= proj_v / proj_component * (sin_val * max_z_error);
        }
        loop.row(i) = bbox_center + proj_v + side_component * side_dir / sin_val;
    }
}

VectorF FastWindingNumberEngine::run(const MatrixFr& queries) {
    const auto num_vertices = m_vertices.rows();
    const auto num_faces = m_faces.rows();
    const auto num_queries = queries.rows();
    using Vector = HDK_Sample::UT_Vector3T<float>;
    using Engine = HDK_Sample::UT_SolidAngle<float, float>;

    std::vector<Vector> vertices(num_vertices);
    for (size_t i=0; i<num_vertices; i++) {
        vertices[i][0] = static_cast<float>(m_vertices(i, 0));
        vertices[i][1] = static_cast<float>(m_vertices(i, 1));
        vertices[i][2] = static_cast<float>(m_vertices(i, 2));
    }

    Engine engine;
    engine.init(num_faces, m_faces.data(), num_vertices, vertices.data());

    VectorF winding_numbers(num_queries);
    for (size_t i=0; i<num_queries; i++) {
        Vector q;
        q[0] = static_cast<float>(queries(i, 0));
        q[1] = static_cast<float>(queries(i, 1));
        q[2] = static_cast<float>(queries(i, 2));
        winding_numbers[i] = engine.computeSolidAngle(q);
    }
    winding_numbers /= 4 * M_PI;

    return winding_numbers;
}

SelfIntersection::SelfIntersection(
        const MatrixFr& vertices, const MatrixIr& faces)
: m_faces(faces) {
    const size_t num_vertices = vertices.rows();
    const size_t dim = vertices.cols();
    const size_t num_faces = faces.rows();
    const size_t vertex_per_face = faces.cols();

    if (dim != 3) {
        throw NotImplementedError(
                "Self intersection check only support 3D");
    }
    if (vertex_per_face != 3) {
        throw NotImplementedError(
                "Self intersection check only works with triangles");
    }

    m_points.resize(num_vertices);
    for (size_t i=0; i<num_vertices; i++) {
        m_points[i] = Point_3(
                vertices(i,0),
                vertices(i,1),
                vertices(i,2));
    }
}

void GeometryCorrectionTable::apply_correction_to_in_plane_edge(
        const Vector3F& edge_dir, MatrixFr& loop) {
    VectorF bbox_min = loop.colwise().minCoeff();
    VectorF bbox_max = loop.colwise().maxCoeff();
    VectorF bbox_center = 0.5 * (bbox_min + bbox_max);

    size_t num_vts = loop.rows();
    size_t dim = loop.cols();
    assert(dim == 3);
    MatrixFr proj_loop(num_vts, dim);

    for (size_t i=0; i<num_vts; i++) {
        const VectorF& v = loop.row(i) - bbox_center.transpose();
        proj_loop.row(i) = Vector3F(v[0], v[1], 0.0);
    }
    Float target_half_height = 1e3; // Something huge to represent inf
    Float target_half_width = proj_loop.row(0).norm();

    Vector2F correction_1 = lookup(target_half_width, target_half_height);
    Vector2F correction_2 = lookup(target_half_height, target_half_width);
    Float half_width  = 0.5 * (correction_1[0] + correction_2[1])
        + 0.05 * num_offset_pixel;
    half_width = std::max(half_width, min_thickness);

    for (size_t i=0; i<num_vts; i++) {
        loop.row(i) += proj_loop.row(i) *
            (-target_half_width + half_width) / target_half_width;
    }
}

GeoMeshPtr GeogramMeshUtils::wire_network_to_geomesh(
        const MatrixFr& vertices, const Matrix2Ir& edges) {
    const size_t dim = vertices.cols();
    const size_t num_vertices = vertices.rows();
    const size_t num_edges = edges.rows();

    auto geo_mesh = std::make_shared<GeoMesh>(dim, false);
    geo_mesh->vertices.clear();
    geo_mesh->vertices.create_vertices(num_vertices);
    geo_mesh->edges.clear();
    geo_mesh->edges.create_edges(num_edges);

    for (size_t i=0; i<num_vertices; i++) {
        auto& p = geo_mesh->vertices.point(i);
        for (size_t j=0; j<dim; j++) {
            p[j] = vertices(i,j);
        }
    }

    for (size_t i=0; i<num_edges; i++) {
        geo_mesh->edges.set_vertex(i, 0, edges(i,0));
        geo_mesh->edges.set_vertex(i, 1, edges(i,1));
    }

    return geo_mesh;
}

GeoMeshPtr GeogramMeshUtils::raw_to_geomesh(
        const MatrixFr& vertices, const MatrixIr& faces) {
    const size_t dim = vertices.cols();
    const size_t vertex_per_face = faces.cols();
    const size_t num_vertices = vertices.rows();
    const size_t num_faces = faces.rows();

    if (vertex_per_face != 3) {
        throw NotImplementedError("Converting non-triangle mesh to "
                "Geogram mesh is not yet implemented");
    }

    auto geo_mesh = std::make_shared<GeoMesh>(dim, false);
    geo_mesh->vertices.clear();
    geo_mesh->vertices.create_vertices(num_vertices);
    geo_mesh->facets.clear();
    geo_mesh->facets.create_triangles(num_faces);

    for (size_t i=0; i<num_vertices; i++) {
        auto& p = geo_mesh->vertices.point(i);
        for (size_t j=0; j<dim; j++) {
            p[j] = vertices(i,j);
        }
    }

    for (size_t i=0; i<num_faces; i++) {
        geo_mesh->facets.set_vertex(i, 0, faces(i,0));
        geo_mesh->facets.set_vertex(i, 1, faces(i,1));
        geo_mesh->facets.set_vertex(i, 2, faces(i,2));
    }

    return geo_mesh;
}

    CarveMeshPtr create_mesh(const MatrixFr& vertices, const MatrixIr& faces) {
        const size_t num_vertices = vertices.rows();
        const size_t num_faces = faces.rows();

        if (vertices.cols() != 3) {
            throw NotImplementedError("Only 3D mesh is supported.");
        }
        if (faces.cols() != 3) {
            throw NotImplementedError("Only triangle mesh is supported.");
        }

        std::vector<CarveVector> points;

        for (size_t i=0; i<num_vertices; i++) {
            const auto& v = vertices.row(i);
            CarveVector p;
            p.v[0] = v[0];
            p.v[1] = v[1];
            p.v[2] = v[2];
            points.push_back(p);
        }

        std::vector<int> raw_faces;
        raw_faces.reserve(num_faces * 4);
        for (size_t i=0; i<num_faces; i++) {
            raw_faces.push_back(3);
            raw_faces.push_back(faces(i,0));
            raw_faces.push_back(faces(i,1));
            raw_faces.push_back(faces(i,2));
        }

        return CarveMeshPtr(new CarveMesh(points, num_faces, raw_faces));
    }

 Polyhedron generate_polyhedron(
         const MatrixFr& vertices, const MatrixIr& faces) {
     Polyhedron P;
     PolyhedronBuilder<HalfedgeDS> triangle(vertices, faces);
     P.delegate(triangle);
     assert(vertices.rows() == P.size_of_vertices());
     assert(faces.rows() == P.size_of_facets());
     return P;
 }

 HashGrid::Ptr compute_vertex_grid(const MatrixFr& vertices, Float cell_size) {
     const size_t dim = vertices.cols();
     const size_t num_vertices = vertices.rows();
     HashGrid::Ptr grid = HashGrid::create(cell_size, dim);
     for (size_t i=0; i<num_vertices; i++) {
         const VectorF& v = vertices.row(i);
         grid->insert(i, v);
     }
     return grid;
 }

void InflatorEngine::save_mesh(const std::string& filename,
        const MatrixFr& vertices, const MatrixIr& faces, VectorF debug) {
    VectorF flattened_vertices(vertices.rows() * vertices.cols());
    std::copy(vertices.data(), vertices.data() + vertices.rows() *
            vertices.cols(), flattened_vertices.data());
    VectorI flattened_faces(faces.rows() * faces.cols());
    std::copy(faces.data(), faces.data() + faces.rows() * faces.cols(),
            flattened_faces.data());
    VectorI voxels = VectorI::Zero(0);

    Mesh::Ptr mesh = MeshFactory().load_data(
            flattened_vertices, flattened_faces, voxels,
            vertices.cols(), faces.cols(), 0).create_shared();
    mesh->add_attribute("debug");
    mesh->set_attribute("debug", debug);

    MeshWriter::Ptr writer = MeshWriter::create(filename);
    writer->with_attribute("debug");
    writer->write_mesh(*mesh);
}

void CGALConvexHull3D::run(const MatrixFr& points) {
    std::list<Point_3> cgal_pts;
    const size_t num_pts = points.rows();
    const size_t dim = points.cols();
    if (dim != 3) {
        std::stringstream err_msg;
        err_msg << "Invalid dim: " << dim << "  Expect dim=3.";
        throw RuntimeError(err_msg.str());
    }

    for (size_t i=0; i<num_pts; i++) {
        const VectorF& p = points.row(i);
        cgal_pts.push_back(Point_3(p[0], p[1], p[2]));
    }

    Polyhedron_3 hull;
    CGAL::convex_hull_3(cgal_pts.begin(), cgal_pts.end(), hull);
    assert(hull.is_closed());
    assert(hull.is_pure_triangle());

    const size_t num_vertices = hull.size_of_vertices();
    const size_t num_faces = hull.size_of_facets();

    m_vertices.resize(num_vertices, dim);
    m_faces.resize(num_faces, 3);

    size_t vertex_count=0;
    for (auto itr=hull.vertices_begin(); itr!=hull.vertices_end(); itr++) {
        const Point_3& p = itr->point();
        m_vertices.coeffRef(vertex_count, 0) = p.x();
        m_vertices.coeffRef(vertex_count, 1) = p.y();
        m_vertices.coeffRef(vertex_count, 2) = p.z();
        itr->id() = vertex_count;
        vertex_count++;
    }

    size_t face_count=0;
    for (auto f_itr=hull.facets_begin(); f_itr!=hull.facets_end(); f_itr++) {
        size_t edge_count=0;
        auto h_itr = f_itr->facet_begin();
        do {
            m_faces.coeffRef(face_count, edge_count) = h_itr->vertex()->id();
            edge_count++;
            h_itr++;
        } while (h_itr != f_itr->facet_begin());
        face_count++;
    }

    compute_index_map(points);
    reorient_faces();
}

    void triangulate(MatrixFr vertices, MatrixIr edges,
            MatrixFr& output_vertices, MatrixIr& output_faces, Float max_area) {
        assert(edges.rows() >= 3);
        MeshCleaner cleaner;
        cleaner.remove_isolated_vertices(vertices, edges);
        cleaner.remove_duplicated_vertices(vertices, edges, 1e-12);
        assert(vertices.rows() >= 3);

        TriangleWrapper triangle(vertices, edges);
        triangle.run(max_area, false, true, true);

        output_vertices = triangle.get_vertices();
        output_faces = triangle.get_faces();
    }

Boundary::Ptr Boundary::extract_surface_boundary_raw(
        MatrixFr& vertices, MatrixIr& faces) {
    VectorF flattened_vertices = Eigen::Map<VectorF>(vertices.data(),
            vertices.rows() * vertices.cols());
    VectorI flattened_faces = Eigen::Map<VectorI>(faces.data(),
            faces.rows() * faces.cols());
    VectorI voxels = VectorI::Zero(0);

    MeshFactory factory;
    Mesh::Ptr mesh = factory.load_data(flattened_vertices, flattened_faces,
            voxels, vertices.cols(), faces.cols(), 0).create();

    return extract_surface_boundary(*mesh);
}

WireNetwork::Ptr MixedMeshTiler::tile() {
    const size_t num_cells = get_num_cells();
    auto transforms = get_tiling_operators();
    auto vars_array = extract_attributes(m_mesh);
    VectorI pattern_id = m_mesh->get_attribute("pattern_id").cast<int>();
    assert(pattern_id.size() == num_cells);

    m_tiled_vertices.clear();
    m_tiled_edges.clear();
    m_tiled_thicknesses.clear();
    m_tiled_offsets.clear();

    size_t v_count = 0;
    auto transform_itr = transforms.begin();
    for (size_t i=0; i<num_cells; i++) {
        set_active_wire_network(pattern_id[i]);
        scale_to_unit_box();
        append_vertices(*transform_itr);
        append_edges(v_count);
        append_thicknesses(vars_array[i]);
        append_offsets(vars_array[i], *transform_itr);

        v_count += m_unit_wire_network->get_num_vertices();
        transform_itr++;
    }

    MatrixFr vertices = vstack(m_tiled_vertices);
    MatrixIr edges = vstack(m_tiled_edges);
    MatrixFr thicknesses = vstack(m_tiled_thicknesses);
    MatrixFr offsets = vstack(m_tiled_offsets);
    assert(edges.minCoeff() >= 0);
    assert(edges.maxCoeff() < vertices.rows());

    WireNetwork::Ptr tiled_network =
        WireNetwork::create_raw(vertices, edges);

    tiled_network->add_attribute("thickness",
            m_target_type == ParameterCommon::VERTEX);
    tiled_network->set_attribute("thickness", thicknesses);
    tiled_network->add_attribute("vertex_offset", true);
    tiled_network->set_attribute("vertex_offset", offsets);

    clean_up(*tiled_network);
    return tiled_network;
}