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

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;
    }
}

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;
    }
}

bool MeshValidation::is_periodic(
        const MatrixFr& vertices, const MatrixIr& faces) {
    const Float EPS = 1e-6;
    HashGrid::Ptr grid = compute_vertex_grid(vertices, EPS);

    Vector3F bbox_min = vertices.colwise().minCoeff();
    Vector3F bbox_max = vertices.colwise().maxCoeff();
    Vector3F bbox_size = bbox_max - bbox_min;

    Vector3F offsets[] = {
        Vector3F( bbox_size[0], 0.0, 0.0),
        Vector3F(-bbox_size[0], 0.0, 0.0),
        Vector3F(0.0, bbox_size[1], 0.0),
        Vector3F(0.0,-bbox_size[1], 0.0),
        Vector3F(0.0, 0.0, bbox_size[2]),
        Vector3F(0.0, 0.0,-bbox_size[2])
    };

    bool result = true;
    const size_t num_vertices = vertices.rows();
    for (size_t i=0; i<num_vertices; i++) {
        const VectorF& v = vertices.row(i);
        if (fabs(v[0] - bbox_min[0]) < EPS) {
            result = result && match(grid, v + offsets[0]);
        }
        if (fabs(v[0] - bbox_max[0]) < EPS) {
            result = result && match(grid, v + offsets[1]);
        }
        if (fabs(v[1] - bbox_min[1]) < EPS) {
            result = result && match(grid, v + offsets[2]);
        }
        if (fabs(v[1] - bbox_max[1]) < EPS) {
            result = result && match(grid, v + offsets[3]);
        }
        if (fabs(v[2] - bbox_min[2]) < EPS) {
            result = result && match(grid, v + offsets[4]);
        }
        if (fabs(v[2] - bbox_max[2]) < EPS) {
            result = result && match(grid, v + offsets[5]);
        }
    }
    return result;
}

VectorI MeshCleaner::compute_importance_level(const MatrixFr& vertices) {
    VectorF bbox_min = vertices.colwise().minCoeff();
    VectorF bbox_max = vertices.colwise().maxCoeff();
    BoxChecker checker(bbox_min, bbox_max);

    const size_t num_vertices = vertices.rows();
    VectorI level = VectorI::Zero(num_vertices);
    for (size_t i=0; i<num_vertices; i++) {
        const VectorF& v = vertices.row(i);
        if (checker.is_on_boundary_corners(v)) {
            level[i] = 3;
        } else if (checker.is_on_boundary_edges(v)) {
            level[i] = 2;
        } else if (checker.is_on_boundary(v)) {
            level[i] = 1;
        }
    }

    return level;
}

void GeometryCorrectionTable::apply_correction_to_out_plane_edge(
        const Vector3F& edge_dir, MatrixFr& loop) {
    const Float EPS = 1e-3;
    assert(fabs(edge_dir[2]) > 0.0);
    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, 3);
    Vector3F proj_edge_dir(edge_dir[0], edge_dir[1], 0.0);

    if (loop.rows() != 4) {
        throw NotImplementedError(
                "Geometry correction supports only square wires");
    }

    for (size_t i=0; i<num_vts; i++) {
        VectorF v = loop.row(i) - bbox_center.transpose();
        Float edge_dir_offset = v[2] / edge_dir[2];
        VectorF proj_v = v - edge_dir * edge_dir_offset;
        proj_loop.row(i) = proj_v;
        assert(fabs(proj_v[2]) < EPS);
    }

    Float dist_01 = (proj_loop.row(0) - proj_loop.row(1)).norm();
    Float dist_12 = (proj_loop.row(1) - proj_loop.row(2)).norm();
    if (dist_01 > dist_12) {
        proj_edge_dir = (proj_loop.row(0) - proj_loop.row(1)) / dist_01;
    } else {
        proj_edge_dir = (proj_loop.row(1) - proj_loop.row(2)) / dist_12;
    }

    const VectorF& corner = proj_loop.row(0);
    Float target_half_height = proj_edge_dir.dot(corner);
    Float target_half_width =
        (corner - proj_edge_dir * target_half_height).norm();
    target_half_height = fabs(target_half_height);

    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;
    Float half_height = 0.5 * (correction_1[1] + correction_2[0])
        + 0.05 * num_offset_pixel;
    half_width = std::max(half_width, min_thickness);
    half_height = std::max(half_height, min_thickness);

    for (size_t i=0; i<num_vts; i++) {
        const VectorF& proj_v = proj_loop.row(i);
        Float height = proj_edge_dir.dot(proj_v);
        VectorF width_dir = (proj_v - proj_edge_dir * height).normalized();
        assert(!isnan(width_dir[0]));
        assert(!isnan(width_dir[1]));
        assert(!isnan(width_dir[2]));

        Float height_sign = (height < 0.0)? -1: 1;
        proj_loop.row(i) = proj_edge_dir * height_sign * half_height
            + width_dir * half_width;
    }

    for (size_t i=0; i<num_vts; i++) {
        const VectorF& proj_v = proj_loop.row(i);
        loop.row(i) = (bbox_center + proj_v - edge_dir *
                edge_dir.dot(proj_v)).transpose();
    }
}

bool MeshValidation::face_source_is_valid(
        const MatrixFr& vertices,
        const MatrixIr& faces,
        const VectorI& face_sources) {
    const Float EPS = 1e-6;
    const size_t num_vertices = vertices.rows();
    const size_t num_faces = faces.rows();

    Vector3F bbox_min = vertices.colwise().minCoeff();
    Vector3F bbox_max = vertices.colwise().maxCoeff();

    bool result = true;
    for (size_t i=0; i<num_faces; i++) {
        const Vector3I& f = faces.row(i);
        const Vector3F& v0 = vertices.row(f[0]);
        const Vector3F& v1 = vertices.row(f[1]);
        const Vector3F& v2 = vertices.row(f[2]);

        if (fabs(v0[0] - bbox_min[0]) < EPS &&
            fabs(v1[0] - bbox_min[0]) < EPS &&
            fabs(v2[0] - bbox_min[0]) < EPS) {
            result = result && (face_sources[i] == 0);
            continue;
        }

        if (fabs(v0[1] - bbox_min[1]) < EPS &&
            fabs(v1[1] - bbox_min[1]) < EPS &&
            fabs(v2[1] - bbox_min[1]) < EPS) {
            result = result && (face_sources[i] == 0);
            continue;
        }

        if (fabs(v0[2] - bbox_min[2]) < EPS &&
            fabs(v1[2] - bbox_min[2]) < EPS &&
            fabs(v2[2] - bbox_min[2]) < EPS) {
            result = result && (face_sources[i] == 0);
            continue;
        }

        if (fabs(v0[0] - bbox_max[0]) < EPS &&
            fabs(v1[0] - bbox_max[0]) < EPS &&
            fabs(v2[0] - bbox_max[0]) < EPS) {
            result = result && (face_sources[i] == 0);
            continue;
        }

        if (fabs(v0[1] - bbox_max[1]) < EPS &&
            fabs(v1[1] - bbox_max[1]) < EPS &&
            fabs(v2[1] - bbox_max[1]) < EPS) {
            result = result && (face_sources[i] == 0);
            continue;
        }

        if (fabs(v0[2] - bbox_max[2]) < EPS &&
            fabs(v1[2] - bbox_max[2]) < EPS &&
            fabs(v2[2] - bbox_max[2]) < EPS) {
            result = result && (face_sources[i] == 0);
            continue;
        }
        result = result && (face_sources[i] != 0);
        if (!result) {
            std::cout << i << ":  ";
            std::cout << face_sources[i] << std::endl;
            std::cout << v0.transpose() << std::endl;
            std::cout << v1.transpose() << std::endl;
            std::cout << v2.transpose() << std::endl;
            return result;
        }
    }
    return result;
}