diff --git a/src/boundary_conditions.cpp b/src/boundary_conditions.cpp
new file mode 100644
index 00000000..40c3eadf
--- /dev/null
+++ b/src/boundary_conditions.cpp
@@ -0,0 +1,67 @@
+#include "default_types.h"
+#include <igl/boundary_conditions.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/eigen/dense.h>
+#include <nanobind/stl/tuple.h>
+
+namespace nb = nanobind;
+using namespace nb::literals;
+namespace pyigl
+{
+  auto boundary_conditions(
+    const nb::DRef<const Eigen::MatrixXN> &V,
+    const nb::DRef<const Eigen::MatrixXI> &Ele,
+    const nb::DRef<const Eigen::MatrixXN> &C,
+    const nb::DRef<const Eigen::VectorXI> &P,
+    const nb::DRef<const Eigen::MatrixXI> &BE,
+    const nb::DRef<const Eigen::MatrixXI> &CE,
+    const nb::DRef<const Eigen::MatrixXI> &CF)
+  {
+
+    Eigen::VectorXI b;
+    Eigen::MatrixXN bc;
+    if(!igl::boundary_conditions(V, Ele, C, P, BE, CE, CF, b, bc))
+    {
+      throw std::runtime_error("boundary_conditions: failed to compute");
+    }
+    return std::make_tuple(b, bc);
+  }
+}
+
+void bind_boundary_conditions(nb::module_ &m)
+{
+  m.def(
+    "boundary_conditions", 
+    &pyigl::boundary_conditions,
+    "V"_a,
+    "Ele"_a,
+    "C"_a,
+    "P"_a=Eigen::VectorXI(),
+    "BE"_a=Eigen::MatrixXI(),
+    "CE"_a=Eigen::MatrixXI(),
+    "CF"_a=Eigen::MatrixXI(),
+R"(Compute boundary conditions for automatic weights computation. This
+function expects that the given mesh (V,Ele) has sufficient samples
+(vertices) exactly at point handle locations and exactly along bone and
+cage edges/faces.
+@param[in] V  #V by dim list of domain vertices
+@param[in] Ele  #Ele by simplex-size list of simplex indices
+@param[in] C  #C by dim list of handle positions
+@param[in] P  #P by 1 list of point handle indices into C
+@param[in] BE  #BE by 2 list of bone edge indices into C
+@param[in] CE  #CE by 2 list of cage edge indices into *P*
+@param[in] CF  #CF by 3 list of (triangular) cage face indices into *P*
+@param[out] b  #b list of boundary indices (indices into V of vertices which have
+    known, fixed values)
+@param[out] bc #b by #weights list of known/fixed values for boundary vertices
+    (notice the #b != #weights in general because #b will include all the
+    intermediary samples along each bone, etc.. The ordering of the
+    weights corresponds to [P;BE]
+throws error if boundary conditions are suspicious:
+  P and BE are empty
+  bc is empty
+  some column of bc doesn't have a 0 (assuming bc has >1 columns)
+  some column of bc doesn't have a 1 (assuming bc has >1 columns))"
+  );
+}
diff --git a/tests/test_all.py b/tests/test_all.py
index f2bab331..9728d16d 100644
--- a/tests/test_all.py
+++ b/tests/test_all.py
@@ -1811,6 +1811,33 @@ def test_resolve_duplicated_faces():
     assert set(map(tuple, F2.tolist())) == {(3, 4, 5)}
 
 
+def test_boundary_conditions():
+    V, F, T = single_tet()
+
+    # Point handles located exactly at every vertex: each handle fixes its own
+    # vertex, so bc is the identity (column j is 1 at handle j, 0 elsewhere).
+    C = V.copy()
+    P = np.arange(V.shape[0], dtype=np.int64)
+    b, bc = igl.boundary_conditions(V, F, C, P)
+    assert b.shape[0] == bc.shape[0]
+    assert bc.shape[1] == P.shape[0]
+    order = np.argsort(b.ravel())
+    assert np.array_equal(b.ravel()[order], P)
+    assert np.allclose(bc[order], np.eye(P.shape[0]))
+
+    # Bone edge with an interior sample: a strip whose middle-bottom vertex (1)
+    # lies exactly on the bone between the two handles at vertices 0 and 2. All
+    # three bottom vertices become boundary samples of the single bone (weight 1).
+    Vb = np.array([[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [1.0, 0.0, 0.0],
+                   [0.0, 1.0, 0.0], [0.5, 1.0, 0.0], [1.0, 1.0, 0.0]])
+    Fb = np.array([[0, 1, 3], [1, 4, 3], [1, 2, 4], [2, 5, 4]], dtype=np.int64)
+    Cb = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
+    BE = np.array([[0, 1]], dtype=np.int64)
+    b, bc = igl.boundary_conditions(Vb, Fb, Cb, BE=BE)
+    assert bc.shape[1] == BE.shape[0]  # one weight column per bone
+    assert set(b.ravel().tolist()) == {0, 1, 2}
+    assert np.allclose(bc, 1.0)
+
 def test_fast_winding_number():
     # Unit-sphere mesh from a subdivided icosahedron.
     V, F = igl.icosahedron()