diff --git a/source/blender/editors/io/io_usd.cc b/source/blender/editors/io/io_usd.cc
index 837791b921a..a0bbcfa6fe2 100644
--- a/source/blender/editors/io/io_usd.cc
+++ b/source/blender/editors/io/io_usd.cc
@@ -429,6 +429,8 @@ static int wm_usd_import_exec(bContext *C, wmOperator *op)
   const bool import_meshes = RNA_boolean_get(op->ptr, "import_meshes");
   const bool import_volumes = RNA_boolean_get(op->ptr, "import_volumes");
   const bool import_shapes = RNA_boolean_get(op->ptr, "import_shapes");
+  const bool import_skeletons = RNA_boolean_get(op->ptr, "import_skeletons");
+  const bool import_blendshapes = RNA_boolean_get(op->ptr, "import_blendshapes");
 
   const bool import_subdiv = RNA_boolean_get(op->ptr, "import_subdiv");
 
@@ -492,6 +494,8 @@ static int wm_usd_import_exec(bContext *C, wmOperator *op)
   params.import_meshes = import_meshes;
   params.import_volumes = import_volumes;
   params.import_shapes = import_shapes;
+  params.import_skeletons = import_skeletons;
+  params.import_blendshapes = import_blendshapes;
   params.prim_path_mask = prim_path_mask;
   params.import_subdiv = import_subdiv;
   params.import_instance_proxies = import_instance_proxies;
@@ -538,6 +542,8 @@ static void wm_usd_import_draw(bContext * /*C*/, wmOperator *op)
   uiItemR(col, ptr, "import_meshes", UI_ITEM_NONE, nullptr, ICON_NONE);
   uiItemR(col, ptr, "import_volumes", UI_ITEM_NONE, nullptr, ICON_NONE);
   uiItemR(col, ptr, "import_shapes", UI_ITEM_NONE, nullptr, ICON_NONE);
+  uiItemR(col, ptr, "import_skeletons", UI_ITEM_NONE, nullptr, ICON_NONE);
+  uiItemR(col, ptr, "import_blendshapes", UI_ITEM_NONE, nullptr, ICON_NONE);
   uiItemR(box, ptr, "prim_path_mask", UI_ITEM_NONE, nullptr, ICON_NONE);
   uiItemR(box, ptr, "scale", UI_ITEM_NONE, nullptr, ICON_NONE);
 
@@ -632,6 +638,8 @@ void WM_OT_usd_import(wmOperatorType *ot)
   RNA_def_boolean(ot->srna, "import_meshes", true, "Meshes", "");
   RNA_def_boolean(ot->srna, "import_volumes", true, "Volumes", "");
   RNA_def_boolean(ot->srna, "import_shapes", true, "Shapes", "");
+  RNA_def_boolean(ot->srna, "import_skeletons", true, "Skeletons", "");
+  RNA_def_boolean(ot->srna, "import_blendshapes", true, "Blend Shapes", "");
 
   RNA_def_boolean(ot->srna,
                   "import_subdiv",
diff --git a/source/blender/io/usd/CMakeLists.txt b/source/blender/io/usd/CMakeLists.txt
index 5880da289a1..a4c867bac8c 100644
--- a/source/blender/io/usd/CMakeLists.txt
+++ b/source/blender/io/usd/CMakeLists.txt
@@ -110,6 +110,8 @@ set(SRC
   intern/usd_reader_stage.cc
   intern/usd_reader_volume.cc
   intern/usd_reader_xform.cc
+  intern/usd_reader_skeleton.cc
+  intern/usd_skel_convert.cc
 
   usd.hh
   usd.h
@@ -142,6 +144,8 @@ set(SRC
   intern/usd_reader_stage.h
   intern/usd_reader_volume.h
   intern/usd_reader_xform.h
+  intern/usd_reader_skeleton.h
+  intern/usd_skel_convert.h
 )
 
 if(WITH_HYDRA)
diff --git a/source/blender/io/usd/intern/usd_capi_import.cc b/source/blender/io/usd/intern/usd_capi_import.cc
index d4ced6009a9..5e98c6b2365 100644
--- a/source/blender/io/usd/intern/usd_capi_import.cc
+++ b/source/blender/io/usd/intern/usd_capi_import.cc
@@ -303,6 +303,10 @@ static void import_startjob(void *customdata, bool *stop, bool *do_update, float
     }
   }
 
+  if (data->params.import_skeletons) {
+    archive->process_armature_modifiers();
+  }
+
   data->import_ok = !data->was_canceled;
 
   *progress = 1.0f;
diff --git a/source/blender/io/usd/intern/usd_reader_mesh.cc b/source/blender/io/usd/intern/usd_reader_mesh.cc
index c308537f508..e9dffe540d3 100644
--- a/source/blender/io/usd/intern/usd_reader_mesh.cc
+++ b/source/blender/io/usd/intern/usd_reader_mesh.cc
@@ -7,6 +7,7 @@
 
 #include "usd_reader_mesh.h"
 #include "usd_reader_material.h"
+#include "usd_skel_convert.h"
 
 #include "BKE_attribute.hh"
 #include "BKE_customdata.h"
@@ -43,6 +44,7 @@
 #include <pxr/usd/usdGeom/primvarsAPI.h>
 #include <pxr/usd/usdGeom/subset.h>
 #include <pxr/usd/usdShade/materialBindingAPI.h>
+#include <pxr/usd/usdSkel/bindingAPI.h>
 
 #include <iostream>
 
@@ -268,6 +270,14 @@ void USDMeshReader::read_object_data(Main *bmain, const double motionSampleTime)
     }
   }
 
+  if (import_params_.import_blendshapes) {
+    import_blendshapes(bmain, object_, prim_);
+  }
+
+  if (import_params_.import_skeletons) {
+    import_mesh_skel_bindings(bmain, object_, prim_);
+  }
+
   USDXformReader::read_object_data(bmain, motionSampleTime);
 }  // namespace blender::io::usd
 
@@ -1099,4 +1109,57 @@ Mesh *USDMeshReader::read_mesh(Mesh *existing_mesh,
   return active_mesh;
 }
 
+std::string USDMeshReader::get_skeleton_path() const
+{
+  /* Make sure we can apply UsdSkelBindingAPI to the prim.
+   * Attempting to apply the API to instance proxies generates
+   * a USD error. */
+  if (!prim_ || prim_.IsInstanceProxy()) {
+    return "";
+  }
+
+  pxr::UsdSkelBindingAPI skel_api = pxr::UsdSkelBindingAPI::Apply(prim_);
+
+  if (!skel_api) {
+    return "";
+  }
+
+  if (pxr::UsdSkelSkeleton skel = skel_api.GetInheritedSkeleton()) {
+    return skel.GetPath().GetAsString();
+  }
+
+  return "";
+}
+
+std::optional<XformResult> USDMeshReader::get_local_usd_xform(const float time) const
+{
+  if (!import_params_.import_skeletons || prim_.IsInstanceProxy()) {
+    /* Use the standard transform computation, since we are ignoring
+     * skinning data. Note that applying the UsdSkelBinding API to an
+     * instance proxy generates a USD error. */
+    return USDXformReader::get_local_usd_xform(time);
+  }
+
+  if (pxr::UsdSkelBindingAPI skel_api = pxr::UsdSkelBindingAPI::Apply(prim_)) {
+    if (skel_api.GetGeomBindTransformAttr().HasAuthoredValue()) {
+      pxr::GfMatrix4d bind_xf;
+      if (skel_api.GetGeomBindTransformAttr().Get(&bind_xf)) {
+        /* The USD bind transform is a matrix of doubles,
+         * but we cast it to GfMatrix4f because Blender expects
+         * a matrix of floats. Also, we assume the transform
+         * is constant over time. */
+        return XformResult(pxr::GfMatrix4f(bind_xf), true);
+      }
+      else {
+        WM_reportf(RPT_WARNING,
+                   "%s: Couldn't compute geom bind transform for %s",
+                   __func__,
+                   prim_.GetPath().GetAsString().c_str());
+      }
+    }
+  }
+
+  return USDXformReader::get_local_usd_xform(time);
+}
+
 }  // namespace blender::io::usd
diff --git a/source/blender/io/usd/intern/usd_reader_mesh.h b/source/blender/io/usd/intern/usd_reader_mesh.h
index e9ce3616c1a..9355c0902d3 100644
--- a/source/blender/io/usd/intern/usd_reader_mesh.h
+++ b/source/blender/io/usd/intern/usd_reader_mesh.h
@@ -56,6 +56,16 @@ class USDMeshReader : public USDGeomReader {
 
   bool topology_changed(const Mesh *existing_mesh, double motionSampleTime) override;
 
+  /**
+   * If the USD mesh prim has a valid UsdSkel schema defined, return the USD path
+   * string to the bound skeleton, if any.  Returns the empty string if no skeleton
+   * binding is defined.
+   *
+   * The returned path is currently used to match armature modifiers with armature
+   * objects during import.
+   */
+  std::string get_skeleton_path() const;
+
  private:
   void process_normals_vertex_varying(Mesh *mesh);
   void process_normals_face_varying(Mesh *mesh);
@@ -82,6 +92,7 @@ class USDMeshReader : public USDGeomReader {
   void read_color_data_primvar(Mesh *mesh,
                                const pxr::UsdGeomPrimvar &color_primvar,
                                const double motionSampleTime);
+
   void read_uv_data_primvar(Mesh *mesh,
                             const pxr::UsdGeomPrimvar &primvar,
                             const double motionSampleTime);
@@ -94,6 +105,12 @@ class USDMeshReader : public USDGeomReader {
                                             const pxr::UsdGeomPrimvar &primvar,
                                             const double motionSampleTime,
                                             MutableSpan<BlenderT> attribute);
+
+  /**
+   * Override transform computation to account for the binding
+   * transformation for skinned meshes.
+   */
+  std::optional<XformResult> get_local_usd_xform(float time) const override;
 };
 
 }  // namespace blender::io::usd
diff --git a/source/blender/io/usd/intern/usd_reader_skeleton.cc b/source/blender/io/usd/intern/usd_reader_skeleton.cc
new file mode 100644
index 00000000000..0fa74e8effd
--- /dev/null
+++ b/source/blender/io/usd/intern/usd_reader_skeleton.cc
@@ -0,0 +1,47 @@
+/* SPDX-FileCopyrightText: 2021 NVIDIA Corporation. All rights reserved.
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "usd_reader_skeleton.h"
+#include "usd_skel_convert.h"
+
+#include "BKE_armature.h"
+#include "BKE_idprop.h"
+#include "BKE_object.h"
+
+#include "DNA_armature_types.h"
+#include "DNA_object_types.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "WM_api.hh"
+
+#include <iostream>
+
+namespace blender::io::usd {
+
+bool USDSkeletonReader::valid() const
+{
+  return skel_ && USDXformReader::valid();
+}
+
+void USDSkeletonReader::create_object(Main *bmain, const double /* motionSampleTime */)
+{
+  object_ = BKE_object_add_only_object(bmain, OB_ARMATURE, name_.c_str());
+
+  bArmature *arm = BKE_armature_add(bmain, name_.c_str());
+  object_->data = arm;
+}
+
+void USDSkeletonReader::read_object_data(Main *bmain, const double motionSampleTime)
+{
+  if (!object_ || !object_->data || !skel_) {
+    return;
+  }
+
+  import_skeleton(bmain, object_, skel_);
+
+  USDXformReader::read_object_data(bmain, motionSampleTime);
+}
+
+}  // namespace blender::io::usd
diff --git a/source/blender/io/usd/intern/usd_reader_skeleton.h b/source/blender/io/usd/intern/usd_reader_skeleton.h
new file mode 100644
index 00000000000..abc6b08e306
--- /dev/null
+++ b/source/blender/io/usd/intern/usd_reader_skeleton.h
@@ -0,0 +1,30 @@
+/* SPDX-FileCopyrightText: 2023 NVIDIA Corporation. All rights reserved.
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+#pragma once
+
+#include "usd.h"
+#include "usd_reader_xform.h"
+
+#include <pxr/usd/usdSkel/skeleton.h>
+
+namespace blender::io::usd {
+
+class USDSkeletonReader : public USDXformReader {
+ private:
+  pxr::UsdSkelSkeleton skel_;
+
+ public:
+  USDSkeletonReader(const pxr::UsdPrim &prim,
+                    const USDImportParams &import_params,
+                    const ImportSettings &settings)
+      : USDXformReader(prim, import_params, settings), skel_(prim)
+  {
+  }
+
+  bool valid() const override;
+  void create_object(Main *bmain, double motionSampleTime) override;
+  void read_object_data(Main *bmain, double motionSampleTime) override;
+};
+
+}  // namespace blender::io::usd
diff --git a/source/blender/io/usd/intern/usd_reader_stage.cc b/source/blender/io/usd/intern/usd_reader_stage.cc
index 317924005d3..7f91507741a 100644
--- a/source/blender/io/usd/intern/usd_reader_stage.cc
+++ b/source/blender/io/usd/intern/usd_reader_stage.cc
@@ -11,6 +11,7 @@
 #include "usd_reader_nurbs.h"
 #include "usd_reader_prim.h"
 #include "usd_reader_shape.h"
+#include "usd_reader_skeleton.h"
 #include "usd_reader_volume.h"
 #include "usd_reader_xform.h"
 
@@ -42,9 +43,12 @@
 #include "BLI_string.h"
 
 #include "BKE_lib_id.h"
+#include "BKE_modifier.h"
 
 #include "DNA_material_types.h"
 
+#include "WM_api.hh"
+
 namespace blender::io::usd {
 
 USDStageReader::USDStageReader(pxr::UsdStageRefPtr stage,
@@ -100,6 +104,9 @@ USDPrimReader *USDStageReader::create_reader_if_allowed(const pxr::UsdPrim &prim
   if (params_.import_volumes && prim.IsA<pxr::UsdVolVolume>()) {
     return new USDVolumeReader(prim, params_, settings_);
   }
+  if (params_.import_skeletons && prim.IsA<pxr::UsdSkelSkeleton>()) {
+    return new USDSkeletonReader(prim, params_, settings_);
+  }
   if (prim.IsA<pxr::UsdGeomImageable>()) {
     return new USDXformReader(prim, params_, settings_);
   }
@@ -134,6 +141,9 @@ USDPrimReader *USDStageReader::create_reader(const pxr::UsdPrim &prim)
   if (prim.IsA<pxr::UsdVolVolume>()) {
     return new USDVolumeReader(prim, params_, settings_);
   }
+  if (prim.IsA<pxr::UsdSkelSkeleton>()) {
+    return new USDSkeletonReader(prim, params_, settings_);
+  }
   if (prim.IsA<pxr::UsdGeomImageable>()) {
     return new USDXformReader(prim, params_, settings_);
   }
@@ -167,6 +177,11 @@ bool USDStageReader::include_by_visibility(const pxr::UsdGeomImageable &imageabl
 
 bool USDStageReader::include_by_purpose(const pxr::UsdGeomImageable &imageable) const
 {
+  if (params_.import_skeletons && imageable.GetPrim().IsA<pxr::UsdSkelSkeleton>()) {
+    /* Always include skeletons, if requested by the user, regardless of purpose. */
+    return true;
+  }
+
   if (params_.import_guide && params_.import_proxy && params_.import_render) {
     /* The options allow any purpose, so we trivially include the prim. */
     return true;
@@ -318,6 +333,48 @@ void USDStageReader::collect_readers(Main *bmain)
   collect_readers(bmain, root);
 }
 
+void USDStageReader::process_armature_modifiers() const
+{
+  /* Iteratate over the skeleton readers to create the
+   * armature object map, which maps a USD skeleton prim
+   * path to the corresponding armature object. */
+  std::map<std::string, Object *> usd_path_to_armature;
+  for (const USDPrimReader *reader : readers_) {
+    if (dynamic_cast<const USDSkeletonReader *>(reader) && reader->object()) {
+      usd_path_to_armature.insert(std::make_pair(reader->prim_path(), reader->object()));
+    }
+  }
+
+  /* Iterate over the mesh readers and set armature objects on armature modifiers. */
+  for (const USDPrimReader *reader : readers_) {
+    if (!reader->object()) {
+      continue;
+    }
+    const USDMeshReader *mesh_reader = dynamic_cast<const USDMeshReader *>(reader);
+    if (!mesh_reader) {
+      continue;
+    }
+    /* Check if the mesh object has an armature modifier. */
+    ModifierData *md = BKE_modifiers_findby_type(reader->object(), eModifierType_Armature);
+    if (!md) {
+      continue;
+    }
+
+    ArmatureModifierData *amd = reinterpret_cast<ArmatureModifierData *>(md);
+
+    /* Assign the armature based on the bound USD skeleton path of the skinned mesh. */
+    std::string skel_path = mesh_reader->get_skeleton_path();
+    std::map<std::string, Object *>::const_iterator it = usd_path_to_armature.find(skel_path);
+    if (it == usd_path_to_armature.end()) {
+      WM_reportf(RPT_WARNING,
+                 "%s: Couldn't find armature object corresponding to USD skeleton %s",
+                 __func__,
+                 skel_path.c_str());
+    }
+    amd->object = it->second;
+  }
+}
+
 void USDStageReader::import_all_materials(Main *bmain)
 {
   BLI_assert(valid());
diff --git a/source/blender/io/usd/intern/usd_reader_stage.h b/source/blender/io/usd/intern/usd_reader_stage.h
index 6104ddb9ba1..b0d50f0902a 100644
--- a/source/blender/io/usd/intern/usd_reader_stage.h
+++ b/source/blender/io/usd/intern/usd_reader_stage.h
@@ -45,6 +45,13 @@ class USDStageReader {
 
   void collect_readers(struct Main *bmain);
 
+  /**
+   * Complete setting up the armature modifiers that
+   * were created for skinned meshes by setting the
+   * modifier object on the corresponding modifier.
+   */
+  void process_armature_modifiers() const;
+
   /* Convert every material prim on the stage to a Blender
    * material, including materials not used by any geometry.
    * Note that collect_readers() must be called before calling
diff --git a/source/blender/io/usd/intern/usd_reader_xform.cc b/source/blender/io/usd/intern/usd_reader_xform.cc
index defde57d567..b9cfa03ce2b 100644
--- a/source/blender/io/usd/intern/usd_reader_xform.cc
+++ b/source/blender/io/usd/intern/usd_reader_xform.cc
@@ -67,37 +67,20 @@ void USDXformReader::read_matrix(float r_mat[4][4] /* local matrix */,
                                  const float scale,
                                  bool *r_is_constant)
 {
-  if (r_is_constant) {
-    *r_is_constant = true;
-  }
+  BLI_assert(r_mat);
+  BLI_assert(r_is_constant);
 
+  *r_is_constant = true;
   unit_m4(r_mat);
 
-  pxr::UsdGeomXformable xformable;
+  std::optional<XformResult> xf_result = get_local_usd_xform(time);
 
-  if (use_parent_xform_) {
-    xformable = pxr::UsdGeomXformable(prim_.GetParent());
-  }
-  else {
-    xformable = pxr::UsdGeomXformable(prim_);
-  }
-
-  if (!xformable) {
-    /* This might happen if the prim is a Scope. */
+  if (!xf_result) {
     return;
   }
 
-  if (r_is_constant) {
-    *r_is_constant = !xformable.TransformMightBeTimeVarying();
-  }
-
-  pxr::GfMatrix4d usd_local_xf;
-  bool reset_xform_stack;
-  xformable.GetLocalTransformation(&usd_local_xf, &reset_xform_stack, time);
-
-  /* Convert the result to a float matrix. */
-  pxr::GfMatrix4f mat4f = pxr::GfMatrix4f(usd_local_xf);
-  mat4f.Get(r_mat);
+  std::get<0>(*xf_result).Get(r_mat);
+  *r_is_constant = std::get<1>(*xf_result);
 
   /* Apply global scaling and rotation only to root objects, parenting
    * will propagate it. */
@@ -168,4 +151,28 @@ bool USDXformReader::is_root_xform_prim() const
   return false;
 }
 
+std::optional<XformResult> USDXformReader::get_local_usd_xform(const float time) const
+{
+  pxr::UsdGeomXformable xformable = use_parent_xform_ ? pxr::UsdGeomXformable(prim_.GetParent()) :
+                                                        pxr::UsdGeomXformable(prim_);
+
+  if (!xformable) {
+    /* This might happen if the prim is a Scope. */
+    return std::nullopt;
+  }
+
+  bool is_constant = !xformable.TransformMightBeTimeVarying();
+
+  bool reset_xform_stack;
+  pxr::GfMatrix4d xform;
+  if (!xformable.GetLocalTransformation(&xform, &reset_xform_stack, time)) {
+    return std::nullopt;
+  }
+
+  /* The USD bind transform is a matrix of doubles,
+   * but we cast it to GfMatrix4f because Blender expects
+   * a matrix of floats. */
+  return XformResult(pxr::GfMatrix4f(xform), is_constant);
+}
+
 }  // namespace blender::io::usd
diff --git a/source/blender/io/usd/intern/usd_reader_xform.h b/source/blender/io/usd/intern/usd_reader_xform.h
index ab357474ca6..41b9cee2d99 100644
--- a/source/blender/io/usd/intern/usd_reader_xform.h
+++ b/source/blender/io/usd/intern/usd_reader_xform.h
@@ -12,6 +12,10 @@
 
 namespace blender::io::usd {
 
+/** A transformation matrix and a boolean indicating
+ * whether the matrix is constant over time. */
+using XformResult = std::tuple<pxr::GfMatrix4f, bool>;
+
 class USDXformReader : public USDPrimReader {
  private:
   bool use_parent_xform_;
@@ -50,6 +54,18 @@ class USDXformReader : public USDPrimReader {
  protected:
   /* Returns true if the contained USD prim is the root of a transform hierarchy. */
   bool is_root_xform_prim() const;
+
+  /**
+   * Return the USD prim's local transformation.
+   *
+   * \param time: Time code for evaluating the transform.
+   *
+   * \return: Optional tuple with the following elements:
+   *          - The transform matrix.
+   *          - A boolean flag indicating whether the matrix
+   *            is constant over time.
+   */
+  virtual std::optional<XformResult> get_local_usd_xform(float time) const;
 };
 
 }  // namespace blender::io::usd
diff --git a/source/blender/io/usd/intern/usd_skel_convert.cc b/source/blender/io/usd/intern/usd_skel_convert.cc
new file mode 100644
index 00000000000..5fa3492cb59
--- /dev/null
+++ b/source/blender/io/usd/intern/usd_skel_convert.cc
@@ -0,0 +1,1038 @@
+/* SPDX-FileCopyrightText: 2023 NVIDIA Corporation. All rights reserved.
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "usd_skel_convert.h"
+
+#include "usd.h"
+
+#include <pxr/usd/usdSkel/animation.h>
+#include <pxr/usd/usdSkel/bindingAPI.h>
+#include <pxr/usd/usdSkel/blendShape.h>
+#include <pxr/usd/usdSkel/cache.h>
+#include <pxr/usd/usdSkel/skeletonQuery.h>
+#include <pxr/usd/usdSkel/utils.h>
+
+#include "DNA_anim_types.h"
+#include "DNA_armature_types.h"
+#include "DNA_key_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_meta_types.h"
+#include "DNA_scene_types.h"
+
+#include "BKE_action.h"
+#include "BKE_armature.h"
+#include "BKE_deform.h"
+#include "BKE_fcurve.h"
+#include "BKE_key.h"
+#include "BKE_lib_id.h"
+#include "BKE_mesh.h"
+#include "BKE_mesh_runtime.hh"
+#include "BKE_modifier.h"
+#include "BKE_object.h"
+#include "BKE_object_deform.h"
+
+#include "BLI_math_vector.h"
+
+#include "ED_armature.hh"
+#include "ED_keyframing.hh"
+#include "ED_mesh.hh"
+
+#include "WM_api.hh"
+
+#include <iostream>
+#include <string>
+#include <vector>
+
+namespace {
+
+/* Utility: return the magnitude of the largest component
+ * of the given vector. */
+inline float max_mag_component(const pxr::GfVec3d &vec)
+{
+  return pxr::GfMax(pxr::GfAbs(vec[0]), pxr::GfAbs(vec[1]), pxr::GfAbs(vec[2]));
+}
+
+/* Utility: create curve at the given array index. */
+FCurve *create_fcurve(const int array_index, const std::string &rna_path)
+{
+  FCurve *fcu = BKE_fcurve_create();
+  fcu->flag = (FCURVE_VISIBLE | FCURVE_SELECTED);
+  fcu->rna_path = BLI_strdup(rna_path.c_str());
+  fcu->array_index = array_index;
+  return fcu;
+}
+
+/* Utility: create curve at the given array index and
+ * add it as a channel to a group. */
+FCurve *create_chan_fcurve(bAction *act,
+                           bActionGroup *grp,
+                           const int array_index,
+                           const std::string &rna_path,
+                           const int totvert)
+{
+  FCurve *fcu = create_fcurve(array_index, rna_path);
+  fcu->totvert = totvert;
+  action_groups_add_channel(act, grp, fcu);
+  return fcu;
+}
+
+/* Utility: add curve sample. */
+void add_bezt(FCurve *fcu,
+              const float frame,
+              const float value,
+              const eBezTriple_Interpolation ipo = BEZT_IPO_LIN)
+{
+  BezTriple bez;
+  memset(&bez, 0, sizeof(BezTriple));
+  bez.vec[1][0] = frame;
+  bez.vec[1][1] = value;
+  bez.ipo = ipo; /* use default interpolation mode here... */
+  bez.f1 = bez.f2 = bez.f3 = SELECT;
+  bez.h1 = bez.h2 = HD_AUTO;
+  insert_bezt_fcurve(fcu, &bez, INSERTKEY_NOFLAGS);
+}
+
+/**
+ * Import a USD skeleton animation as an action on the given armature object.
+ * This assumes bones have already been created on the armature.
+ *
+ * \param bmain: Main pointer
+ * \param arm_obj: Armature object to which the action will be added
+ * \param skel_query: The USD skeleton query for reading the animation
+ * \param joint_to_bone_map: Map a USD skeleton joint name to a bone name
+ */
+void import_skeleton_curves(Main *bmain,
+                            Object *arm_obj,
+                            const pxr::UsdSkelSkeletonQuery &skel_query,
+                            const std::map<pxr::TfToken, std::string> &joint_to_bone_map)
+
+{
+  if (!(bmain && arm_obj && skel_query)) {
+    return;
+  }
+
+  if (joint_to_bone_map.empty()) {
+    return;
+  }
+
+  const pxr::UsdSkelAnimQuery &anim_query = skel_query.GetAnimQuery();
+
+  if (!anim_query) {
+    /* No animation is defined. */
+    return;
+  }
+
+  std::vector<double> samples;
+  anim_query.GetJointTransformTimeSamples(&samples);
+
+  if (samples.empty()) {
+    return;
+  }
+
+  const size_t num_samples = samples.size();
+
+  /* Create the action on the armature. */
+  bAction *act = ED_id_action_ensure(bmain, (ID *)&arm_obj->id);
+
+  /* Create the curves. */
+
+  /* Get the joint paths. */
+  pxr::VtTokenArray joint_order = skel_query.GetJointOrder();
+
+  std::vector<FCurve *> loc_curves;
+  std::vector<FCurve *> rot_curves;
+  std::vector<FCurve *> scale_curves;
+
+  /* Iterate over the joints and create the corresponding curves for the bones. */
+  for (const pxr::TfToken &joint : joint_order) {
+    std::map<pxr::TfToken, std::string>::const_iterator it = joint_to_bone_map.find(joint);
+
+    if (it == joint_to_bone_map.end()) {
+      /* This joint doesn't correspond to any bone we created.
+       * Add null placeholders for the channel curves. */
+      loc_curves.push_back(nullptr);
+      loc_curves.push_back(nullptr);
+      loc_curves.push_back(nullptr);
+      rot_curves.push_back(nullptr);
+      rot_curves.push_back(nullptr);
+      rot_curves.push_back(nullptr);
+      rot_curves.push_back(nullptr);
+      scale_curves.push_back(nullptr);
+      scale_curves.push_back(nullptr);
+      scale_curves.push_back(nullptr);
+      continue;
+    }
+
+    bActionGroup *grp = action_groups_add_new(act, it->second.c_str());
+
+    /* Add translation curves. */
+    std::string rna_path = "pose.bones[\"" + it->second + "\"].location";
+    loc_curves.push_back(create_chan_fcurve(act, grp, 0, rna_path, num_samples));
+    loc_curves.push_back(create_chan_fcurve(act, grp, 1, rna_path, num_samples));
+    loc_curves.push_back(create_chan_fcurve(act, grp, 2, rna_path, num_samples));
+
+    /* Add rotation curves. */
+    rna_path = "pose.bones[\"" + it->second + "\"].rotation_quaternion";
+    rot_curves.push_back(create_chan_fcurve(act, grp, 0, rna_path, num_samples));
+    rot_curves.push_back(create_chan_fcurve(act, grp, 1, rna_path, num_samples));
+    rot_curves.push_back(create_chan_fcurve(act, grp, 2, rna_path, num_samples));
+    rot_curves.push_back(create_chan_fcurve(act, grp, 3, rna_path, num_samples));
+
+    /* Add scale curves. */
+    rna_path = "pose.bones[\"" + it->second + "\"].scale";
+    scale_curves.push_back(create_chan_fcurve(act, grp, 0, rna_path, num_samples));
+    scale_curves.push_back(create_chan_fcurve(act, grp, 1, rna_path, num_samples));
+    scale_curves.push_back(create_chan_fcurve(act, grp, 2, rna_path, num_samples));
+  }
+
+  /* Sanity checks: make sure we have a curve entry for each joint. */
+  if (loc_curves.size() != joint_order.size() * 3) {
+    std::cout << "PROGRAMMER ERROR: location curve count mismatch\n";
+    return;
+  }
+
+  if (rot_curves.size() != joint_order.size() * 4) {
+    std::cout << "PROGRAMMER ERROR: rotation curve count mismatch\n";
+    return;
+  }
+
+  if (scale_curves.size() != joint_order.size() * 3) {
+    std::cout << "PROGRAMMER ERROR: scale curve count mismatch\n";
+    return;
+  }
+
+  /* The curve for each joint represents the transform relative
+   * to the bind transform in joint-local space. I.e.,
+   *
+   * jointLocalTransform * inv(jointLocalBindTransform)
+   *
+   * There doesn't appear to be a way to query the joint-local
+   * bind transform through the API, so we have to compute it
+   * ourselves from the world bind transforms and the skeleton
+   * topology.
+   */
+
+  /* Get the world space joint transforms at bind time. */
+  pxr::VtMatrix4dArray bind_xforms;
+  if (!skel_query.GetJointWorldBindTransforms(&bind_xforms)) {
+    WM_reportf(RPT_WARNING,
+               "%s: Couldn't get world bind transforms for skeleton %s",
+               __func__,
+               skel_query.GetSkeleton().GetPrim().GetPath().GetAsString().c_str());
+    return;
+  }
+
+  if (bind_xforms.size() != joint_order.size()) {
+    WM_reportf(RPT_WARNING,
+               "%s: Number of bind transforms doesn't match the number of joints for skeleton %s",
+               __func__,
+               skel_query.GetSkeleton().GetPrim().GetPath().GetAsString().c_str());
+    return;
+  }
+
+  const pxr::UsdSkelTopology &skel_topology = skel_query.GetTopology();
+
+  pxr::VtMatrix4dArray joint_local_bind_xforms(bind_xforms.size());
+  for (int i = 0; i < bind_xforms.size(); ++i) {
+    const int parent_id = skel_topology.GetParent(i);
+
+    if (parent_id >= 0) {
+      /* This is a non-root joint.  Compute the bind transform of the joint
+       * relative to its parent. */
+      joint_local_bind_xforms[i] = bind_xforms[i] * bind_xforms[parent_id].GetInverse();
+    }
+    else {
+      /* This is the root joint. */
+      joint_local_bind_xforms[i] = bind_xforms[i];
+    }
+  }
+
+  /* Set the curve samples. */
+  for (const double frame : samples) {
+    pxr::VtMatrix4dArray joint_local_xforms;
+    if (!skel_query.ComputeJointLocalTransforms(&joint_local_xforms, frame)) {
+      std::cout << "WARNING: couldn't compute joint local transforms on frame " << frame
+                << std::endl;
+      continue;
+    }
+
+    if (joint_local_xforms.size() != joint_order.size()) {
+      std::cout << "WARNING: number of joint local transform entries " << joint_local_xforms.size()
+                << " doesn't match the number of joints " << joint_order.size() << std::endl;
+      continue;
+    }
+
+    for (int i = 0; i < joint_local_xforms.size(); ++i) {
+      pxr::GfMatrix4d bone_xform = joint_local_xforms[i] * joint_local_bind_xforms[i].GetInverse();
+
+      pxr::GfVec3f t;
+      pxr::GfQuatf qrot;
+      pxr::GfVec3h s;
+
+      if (!pxr::UsdSkelDecomposeTransform(bone_xform, &t, &qrot, &s)) {
+        std::cout << "WARNING: error decomposing matrix on frame " << frame << std::endl;
+        continue;
+      }
+
+      const float re = qrot.GetReal();
+      const pxr::GfVec3f &im = qrot.GetImaginary();
+
+      for (int j = 0; j < 3; ++j) {
+        const int k = 3 * i + j;
+        if (k >= loc_curves.size()) {
+          std::cout << "PROGRAMMER ERROR: out of bounds translation curve index." << std::endl;
+          break;
+        }
+        if (FCurve *fcu = loc_curves[k]) {
+          add_bezt(fcu, frame, t[j]);
+        }
+      }
+
+      for (int j = 0; j < 4; ++j) {
+        const int k = 4 * i + j;
+        if (k >= rot_curves.size()) {
+          std::cout << "PROGRAMMER ERROR: out of bounds rotation curve index." << std::endl;
+          break;
+        }
+        if (FCurve *fcu = rot_curves[k]) {
+          if (j == 0) {
+            add_bezt(fcu, frame, re);
+          }
+          else {
+            add_bezt(fcu, frame, im[j - 1]);
+          }
+        }
+      }
+
+      for (int j = 0; j < 3; ++j) {
+        const int k = 3 * i + j;
+        if (k >= scale_curves.size()) {
+          std::cout << "PROGRAMMER ERROR: out of bounds scale curve index." << std::endl;
+          break;
+        }
+        if (FCurve *fcu = scale_curves[k]) {
+          add_bezt(fcu, frame, s[j]);
+        }
+      }
+    }
+  }
+
+  /* Recalculate curve handles. */
+  auto recalc_handles = [](FCurve *fcu) { BKE_fcurve_handles_recalc(fcu); };
+  std::for_each(loc_curves.begin(), loc_curves.end(), recalc_handles);
+  std::for_each(rot_curves.begin(), rot_curves.end(), recalc_handles);
+  std::for_each(scale_curves.begin(), scale_curves.end(), recalc_handles);
+}
+
+}  // End anonymous namespace.
+
+namespace blender::io::usd {
+
+void import_blendshapes(Main *bmain,
+                        Object *mesh_obj,
+                        const pxr::UsdPrim &prim,
+                        const bool import_anim)
+{
+  if (!(mesh_obj && mesh_obj->data && mesh_obj->type == OB_MESH && prim)) {
+    return;
+  }
+
+  if (prim.IsInstanceProxy()) {
+    /* Attempting to create a UsdSkelBindingAPI for
+     * instance proxies generates USD errors. */
+    return;
+  }
+
+  pxr::UsdSkelBindingAPI skel_api = pxr::UsdSkelBindingAPI::Apply(prim);
+
+  if (!skel_api) {
+    /* No skel binding. */
+    return;
+  }
+
+  /* Get the blend shape targets, which are the USD paths to the
+   * blend shape primitives. */
+
+  if (!skel_api.GetBlendShapeTargetsRel().HasAuthoredTargets()) {
+    /* No targets. */
+    return;
+  }
+
+  pxr::SdfPathVector targets;
+  if (!skel_api.GetBlendShapeTargetsRel().GetTargets(&targets)) {
+    WM_reportf(RPT_WARNING,
+               "%s: Couldn't get blendshape targets for prim %s",
+               __func__,
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  if (targets.empty()) {
+    return;
+  }
+
+  if (!skel_api.GetBlendShapesAttr().HasAuthoredValue()) {
+    return;
+  }
+
+  /* Get the blend shape name tokens. */
+  pxr::VtTokenArray blendshapes;
+  if (!skel_api.GetBlendShapesAttr().Get(&blendshapes)) {
+    return;
+  }
+
+  if (blendshapes.empty()) {
+    return;
+  }
+
+  /* Sanity check. */
+  if (targets.size() != blendshapes.size()) {
+    WM_reportf(RPT_WARNING,
+               "%s: Number of blendshapes doesn't match number of blendshape targets for prim %s",
+               __func__,
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  pxr::UsdStageRefPtr stage = prim.GetStage();
+
+  if (!stage) {
+    WM_reportf(RPT_WARNING,
+               "%s: Couldn't get stage for prim %s",
+               __func__,
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  Mesh *mesh = static_cast<Mesh *>(mesh_obj->data);
+
+  /* Insert key to source mesh. */
+  Key *key = BKE_key_add(bmain, (ID *)mesh);
+  key->type = KEY_RELATIVE;
+
+  mesh->key = key;
+
+  /* Insert basis key. */
+  KeyBlock *kb = BKE_keyblock_add(key, "Basis");
+  BKE_keyblock_convert_from_mesh(mesh, key, kb);
+
+  /* Keep track of the shapkeys we're adding, for
+   * validation when creating curves later. */
+  std::set<pxr::TfToken> shapekey_names;
+
+  for (int i = 0; i < targets.size(); ++i) {
+    /* Get USD path to blend shape. */
+    const pxr::SdfPath &path = targets[i];
+    pxr::UsdSkelBlendShape blendshape(stage->GetPrimAtPath(path));
+
+    if (!blendshape) {
+      continue;
+    }
+
+    /* Get the blend shape offests. */
+    if (!blendshape.GetOffsetsAttr().HasAuthoredValue()) {
+      /* Blend shape has no authored offsets. */
+      continue;
+    }
+
+    pxr::VtVec3fArray offsets;
+    if (!blendshape.GetOffsetsAttr().Get(&offsets)) {
+      WM_reportf(RPT_WARNING,
+                 "%s: Couldn't get offsets for blend shape %s",
+                 __func__,
+                 path.GetAsString().c_str());
+      continue;
+    }
+
+    if (offsets.empty()) {
+      WM_reportf(
+          RPT_WARNING, "%s: No offsets for blend shape %s", __func__, path.GetAsString().c_str());
+      continue;
+    }
+
+    shapekey_names.insert(blendshapes[i]);
+
+    /* Add the key block. */
+    kb = BKE_keyblock_add(key, blendshapes[i].GetString().c_str());
+    BKE_keyblock_convert_from_mesh(mesh, key, kb);
+
+    /* if authored, point indices are indices into the original mesh
+     * that correspond to the values in the offsets array. */
+    pxr::VtArray<int> point_indices;
+    if (blendshape.GetPointIndicesAttr().HasAuthoredValue()) {
+      blendshape.GetPointIndicesAttr().Get(&point_indices);
+    }
+
+    float *fp = static_cast<float *>(kb->data);
+
+    if (point_indices.empty()) {
+      /* Iterate over all key block elements and add the corresponding
+       * offset to the key block point. */
+      for (int a = 0; a < kb->totelem; ++a, fp += 3) {
+        if (a >= offsets.size()) {
+          WM_reportf(
+              RPT_WARNING,
+              "%s: Number of offsets greater than number of mesh vertices for blend shape %s",
+              __func__,
+              path.GetAsString().c_str());
+          break;
+        }
+        add_v3_v3(fp, offsets[a].data());
+      }
+    }
+    else {
+      /* Iterate over the point indices and add the offset to the corresponding
+       * key block point. */
+      int a = 0;
+      for (int i : point_indices) {
+        if (i < 0 || i > kb->totelem) {
+          std::cerr << "Out of bounds point index " << i << " for blendshape " << path
+                    << std::endl;
+          ++a;
+          continue;
+        }
+        if (a >= offsets.size()) {
+          WM_reportf(
+              RPT_WARNING,
+              "%s: Number of offsets greater than number of mesh vertices for blend shape %s",
+              __func__,
+              path.GetAsString().c_str());
+          break;
+        }
+        add_v3_v3(&fp[3 * i], offsets[a].data());
+        ++a;
+      }
+    }
+  }
+
+  if (!import_anim) {
+    /* We're not importing animation, so we are done. */
+    return;
+  }
+
+  /* Get the blend animation source from the skeleton. */
+
+  pxr::UsdSkelSkeleton skel_prim = skel_api.GetInheritedSkeleton();
+
+  if (!skel_prim) {
+    return;
+  }
+
+  skel_api = pxr::UsdSkelBindingAPI::Apply(skel_prim.GetPrim());
+
+  if (!skel_api) {
+    return;
+  }
+
+  pxr::UsdPrim anim_prim = skel_api.GetInheritedAnimationSource();
+
+  if (!anim_prim) {
+    return;
+  }
+
+  pxr::UsdSkelAnimation skel_anim(anim_prim);
+
+  if (!skel_anim) {
+    return;
+  }
+
+  /* Check if a blend shape weight animation was authored. */
+  if (!skel_anim.GetBlendShapesAttr().HasAuthoredValue()) {
+    return;
+  }
+
+  pxr::UsdAttribute weights_attr = skel_anim.GetBlendShapeWeightsAttr();
+
+  if (!(weights_attr && weights_attr.HasAuthoredValue())) {
+    return;
+  }
+
+  /* Get the animation time samples. */
+  std::vector<double> times;
+  if (!weights_attr.GetTimeSamples(&times)) {
+    return;
+  }
+
+  if (times.empty()) {
+    return;
+  }
+
+  /* Get the blend shape name tokens. */
+  if (!skel_anim.GetBlendShapesAttr().Get(&blendshapes)) {
+    return;
+  }
+
+  if (blendshapes.empty()) {
+    return;
+  }
+
+  const size_t num_samples = times.size();
+
+  /* Create the animation and curves. */
+  bAction *act = ED_id_action_ensure(bmain, (ID *)&key->id);
+  std::vector<FCurve *> curves;
+
+  for (auto blendshape_name : blendshapes) {
+    if (shapekey_names.find(blendshape_name) == shapekey_names.end()) {
+      /* We didn't create a shapekey fo this blendshape, so we don't
+       * create a curve and insert a null placeholder in the curve array. */
+      curves.push_back(nullptr);
+      continue;
+    }
+
+    /* Create the curve for this shape key. */
+    std::string rna_path = "key_blocks[\"" + blendshape_name.GetString() + "\"].value";
+    FCurve *fcu = create_fcurve(0, rna_path);
+    fcu->totvert = num_samples;
+    curves.push_back(fcu);
+    BLI_addtail(&act->curves, fcu);
+  }
+
+  /* Add the weight time samples to the curves. */
+  for (double frame : times) {
+    pxr::VtFloatArray weights;
+    if (!weights_attr.Get(&weights, frame)) {
+      std::cerr << "Couldn't get blendshape weights for time " << frame << std::endl;
+      continue;
+    }
+
+    if (weights.size() != curves.size()) {
+      std::cerr << "Programmer error: number of weight samples doesn't match number of shapekey "
+                   "curve entries for frame "
+                << frame << std::endl;
+      continue;
+    }
+
+    for (int wi = 0; wi < weights.size(); ++wi) {
+      if (curves[wi] != nullptr) {
+        add_bezt(curves[wi], frame, weights[wi]);
+      }
+    }
+  }
+
+  /* Recalculate curve handles. */
+  auto recalc_handles = [](FCurve *fcu) { BKE_fcurve_handles_recalc(fcu); };
+  std::for_each(curves.begin(), curves.end(), recalc_handles);
+}
+
+void import_skeleton(Main *bmain,
+                     Object *arm_obj,
+                     const pxr::UsdSkelSkeleton &skel,
+                     const bool import_anim)
+{
+  if (!(arm_obj && arm_obj->data && arm_obj->type == OB_ARMATURE)) {
+    return;
+  }
+
+  pxr::UsdSkelCache skel_cache;
+  pxr::UsdSkelSkeletonQuery skel_query = skel_cache.GetSkelQuery(skel);
+
+  if (!skel_query.IsValid()) {
+    WM_reportf(RPT_WARNING,
+               "%s: Couldn't query skeleton %s",
+               __func__,
+               skel.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  const pxr::UsdSkelTopology &skel_topology = skel_query.GetTopology();
+
+  pxr::VtTokenArray joint_order = skel_query.GetJointOrder();
+
+  if (joint_order.size() != skel_topology.size()) {
+    WM_reportf(RPT_WARNING,
+               "%s: Topology and joint order size mismatch for skeleton %s",
+               __func__,
+               skel.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  bArmature *arm = static_cast<bArmature *>(arm_obj->data);
+
+  /* Set the armature to edit mode when creating the bones. */
+  ED_armature_to_edit(arm);
+
+  /* The bones we create, stored in the skeleton's joint order. */
+  std::vector<EditBone *> edit_bones;
+
+  /* Keep track of the bones we create for each joint.
+   * We'll need this when creating animation curves
+   * later. */
+  std::map<pxr::TfToken, std::string> joint_to_bone_map;
+
+  /* Create the bones. */
+  for (const pxr::TfToken &joint : joint_order) {
+    std::string name = pxr::SdfPath(joint).GetName();
+    EditBone *bone = ED_armature_ebone_add(arm, name.c_str());
+    if (!bone) {
+      WM_reportf(
+          RPT_WARNING, "%s: Couldn't add bone for joint %s", __func__, joint.GetString().c_str());
+      edit_bones.push_back(nullptr);
+      continue;
+    }
+    joint_to_bone_map.insert(std::make_pair(joint, bone->name));
+    edit_bones.push_back(bone);
+  }
+
+  /* Sanity check: we should have created a bone for each joint. */
+  const size_t num_joints = skel_topology.GetNumJoints();
+  if (edit_bones.size() != num_joints) {
+    WM_reportf(RPT_WARNING,
+               "%s: Mismatch in bone and joint counts for skeleton %s",
+               __func__,
+               skel.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  /* Get the world space joint transforms at bind time. */
+  pxr::VtMatrix4dArray bind_xforms;
+  if (!skel_query.GetJointWorldBindTransforms(&bind_xforms)) {
+    WM_reportf(RPT_WARNING,
+               "%s: Couldn't get world bind transforms for skeleton %s",
+               __func__,
+               skel.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  if (bind_xforms.size() != num_joints) {
+    WM_reportf(RPT_WARNING,
+               "%s:  Mismatch in bind xforms and joint counts for skeleton %s",
+               __func__,
+               skel.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  /* Check if any bone matrices have negative determinants,
+   * indicating negative scales, possibly due to mirroring
+   * operations.  Such matrices can't be propery converted
+   * to Blender's axis/roll bone representation (see
+   * https://projects.blender.org/blender/blender/issues/82930).
+   * If we detect such matrices, we will flag an error and won't
+   * try to import the animation, since the rotations would
+   * be incorrect in such cases.  Unfortunately, the Pixar
+   * UsdSkel examples of the "HumanFemale" suffer from
+   * this issue. */
+  bool negative_determinant = false;
+
+  /* Set bone rest transforms. */
+  for (size_t i = 0; i < num_joints; ++i) {
+    EditBone *ebone = edit_bones[i];
+
+    if (!ebone) {
+      continue;
+    }
+
+    pxr::GfMatrix4f mat(bind_xforms[i]);
+
+    float mat4[4][4];
+    mat.Get(mat4);
+
+    pxr::GfVec3f head(0.0f, 0.0f, 0.0f);
+    pxr::GfVec3f tail(0.0f, 1.0f, 0.0f);
+
+    copy_v3_v3(ebone->head, head.data());
+    copy_v3_v3(ebone->tail, tail.data());
+
+    ED_armature_ebone_from_mat4(ebone, mat4);
+
+    if (mat.GetDeterminant() < 0.0) {
+      negative_determinant = true;
+    }
+  }
+
+  bool valid_skeleton = true;
+  if (negative_determinant) {
+    valid_skeleton = false;
+    WM_reportf(RPT_WARNING,
+               "USD Skeleton Import: bone matrices with negative determinants detected in prim %s."
+               "Such matrices may indicate negative scales, possibly due to mirroring operations, "
+               "and can't currently be converted to Blender's bone representation.  "
+               "The skeletal animation won't be imported",
+               skel.GetPath().GetAsString().c_str());
+  }
+
+  /* Set bone parenting.  In addition, scale bones to account
+   * for separation between parents and children, so that the
+   * bone size is in proportion with the overall skeleton hierarchy.
+   * USD skeletons are composed of joints which we imperfectly
+   * represent as bones. */
+
+  /* This will record the child bone indices per parent bone,
+   * to simplify accessing children when computing lengths. */
+  std::vector<std::vector<int>> child_bones(num_joints);
+
+  for (size_t i = 0; i < num_joints; ++i) {
+    const int parent_idx = skel_topology.GetParent(i);
+    if (parent_idx < 0) {
+      continue;
+    }
+    if (parent_idx >= edit_bones.size()) {
+      std::cout << "WARNING: out of bounds parent index for bone " << pxr::SdfPath(joint_order[i])
+                << " for skeleton " << skel.GetPath() << std::endl;
+      continue;
+    }
+
+    child_bones[parent_idx].push_back(i);
+    if (edit_bones[i] && edit_bones[parent_idx]) {
+      edit_bones[i]->parent = edit_bones[parent_idx];
+    }
+  }
+
+  float avg_len_scale = 0;
+  for (size_t i = 0; i < num_joints; ++i) {
+
+    /* If the bone has any children, scale its length
+     * by the distance between this bone's head
+     * and the average head location of its children. */
+
+    if (child_bones[i].empty()) {
+      continue;
+    }
+
+    EditBone *parent = edit_bones[i];
+    if (!parent) {
+      continue;
+    }
+
+    pxr::GfVec3f avg_child_head(0);
+    for (int j : child_bones[i]) {
+      EditBone *child = edit_bones[j];
+      if (!child) {
+        continue;
+      }
+      pxr::GfVec3f child_head(child->head);
+      avg_child_head += child_head;
+    }
+
+    avg_child_head /= child_bones[i].size();
+
+    pxr::GfVec3f parent_head(parent->head);
+    pxr::GfVec3f parent_tail(parent->tail);
+
+    const float new_len = (avg_child_head - parent_head).GetLength();
+
+    /* Check for epsilon relative to the parent head before scaling. */
+    if (new_len > .00001 * max_mag_component(parent_head)) {
+      parent_tail = parent_head + (parent_tail - parent_head).GetNormalized() * new_len;
+      copy_v3_v3(parent->tail, parent_tail.data());
+      avg_len_scale += new_len;
+    }
+  }
+
+  /* Scale terminal bones by the average length scale. */
+  avg_len_scale /= num_joints;
+
+  for (size_t i = 0; i < num_joints; ++i) {
+    if (!child_bones[i].empty()) {
+      /* Not a terminal bone. */
+      continue;
+    }
+    EditBone *bone = edit_bones[i];
+    if (!bone) {
+      continue;
+    }
+    pxr::GfVec3f head(bone->head);
+
+    /* Check for epsilon relative to the head before scaling. */
+    if (avg_len_scale > .00001 * max_mag_component(head)) {
+      pxr::GfVec3f tail(bone->tail);
+      tail = head + (tail - head).GetNormalized() * avg_len_scale;
+      copy_v3_v3(bone->tail, tail.data());
+    }
+  }
+
+  /* Get out of edit mode. */
+  ED_armature_from_edit(bmain, arm);
+  ED_armature_edit_free(arm);
+
+  if (import_anim && valid_skeleton) {
+    import_skeleton_curves(bmain, arm_obj, skel_query, joint_to_bone_map);
+  }
+}
+
+void import_mesh_skel_bindings(Main *bmain, Object *mesh_obj, const pxr::UsdPrim &prim)
+{
+  if (!(bmain && mesh_obj && mesh_obj->type == OB_MESH && prim)) {
+    return;
+  }
+
+  if (prim.IsInstanceProxy()) {
+    /* Attempting to create a UsdSkelBindingAPI for
+     * instance proxies generates USD errors. */
+    return;
+  }
+
+  pxr::UsdSkelBindingAPI skel_api = pxr::UsdSkelBindingAPI::Apply(prim);
+
+  if (!skel_api) {
+    return;
+  }
+
+  pxr::UsdSkelSkeleton skel = skel_api.GetInheritedSkeleton();
+
+  if (!skel) {
+    return;
+  }
+
+  /* Get the joint identifiers from the skeleton. We will
+   * need these to construct deform groups. */
+  pxr::VtArray<pxr::TfToken> joints;
+
+  if (skel_api.GetJointsAttr().HasAuthoredValue()) {
+    skel_api.GetJointsAttr().Get(&joints);
+  }
+  else if (skel.GetJointsAttr().HasAuthoredValue()) {
+    skel.GetJointsAttr().Get(&joints);
+  }
+
+  if (joints.empty()) {
+    return;
+  }
+
+  /* Get the joint indices, which specify which joints influence a given point. */
+  pxr::UsdGeomPrimvar joint_indices_primvar = skel_api.GetJointIndicesPrimvar();
+  if (!(joint_indices_primvar && joint_indices_primvar.HasAuthoredValue())) {
+    return;
+  }
+
+  /* Get the weights, which specify the weight of a joint on a given point. */
+  pxr::UsdGeomPrimvar joint_weights_primvar = skel_api.GetJointWeightsPrimvar();
+  if (!(joint_weights_primvar && joint_weights_primvar.HasAuthoredValue())) {
+    return;
+  }
+
+  /* Element size specifies the number of joints that might influece a given point.
+   * This is the stride we take when accessing the indices and weights for a
+   * given point. */
+  int joint_indices_elem_size = joint_indices_primvar.GetElementSize();
+  int joint_weights_elem_size = joint_weights_primvar.GetElementSize();
+
+  /* We expect the element counts to match. */
+  if (joint_indices_elem_size != joint_weights_elem_size) {
+    WM_reportf(RPT_WARNING,
+               "%s: Joint weights and joint indices element size mismatch for prim %s",
+               __func__,
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  /* Get the joint indices and weights. */
+  pxr::VtIntArray joint_indices;
+  joint_indices_primvar.ComputeFlattened(&joint_indices);
+
+  pxr::VtFloatArray joint_weights;
+  joint_weights_primvar.ComputeFlattened(&joint_weights);
+
+  if (joint_indices.empty() || joint_weights.empty()) {
+    return;
+  }
+
+  if (joint_indices.size() != joint_weights.size()) {
+    WM_reportf(RPT_WARNING,
+               "%s: Joint weights and joint indices size mismatch size mismatch for prim %s",
+               __func__,
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  Mesh *mesh = static_cast<Mesh *>(mesh_obj->data);
+
+  const pxr::TfToken interp = joint_weights_primvar.GetInterpolation();
+
+  /* Sanity check: we expect only vertex or constant interpolation. */
+  if (interp != pxr::UsdGeomTokens->vertex && interp != pxr::UsdGeomTokens->constant) {
+    WM_reportf(RPT_WARNING,
+               "%s: Unexpected joint weights interpolation type %s for prim %s",
+               __func__,
+               interp.GetString().c_str(),
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  /* Sanity check: make sure we have the expected number of values for the interpolation type. */
+  if (interp == pxr::UsdGeomTokens->vertex &&
+      joint_weights.size() != mesh->totvert * joint_weights_elem_size)
+  {
+    WM_reportf(RPT_WARNING,
+               "%s: Joint weights of unexpected size for vertex interpolation for prim %s",
+               __func__,
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  if (interp == pxr::UsdGeomTokens->constant && joint_weights.size() != joint_weights_elem_size) {
+    WM_reportf(RPT_WARNING,
+               "%s: Joint weights of unexpected size for constant interpolation for prim %s",
+               __func__,
+               prim.GetPath().GetAsString().c_str());
+    return;
+  }
+
+  /* Determine which joint indices are used for skinning this prim. */
+  std::vector<int> used_indices;
+  for (int index : joint_indices) {
+    if (std::find(used_indices.begin(), used_indices.end(), index) == used_indices.end()) {
+      /* We haven't accounted for this index yet. */
+      if (index < 0 || index >= joints.size()) {
+        std::cerr << "Out of bound joint index " << index << std::endl;
+        continue;
+      }
+      used_indices.push_back(index);
+    }
+  }
+
+  if (used_indices.empty()) {
+    return;
+  }
+
+  if (BKE_object_defgroup_data_create(static_cast<ID *>(mesh_obj->data)) == NULL) {
+    WM_reportf(RPT_WARNING,
+               "%s: Error creating deform group data for mesh %s",
+               __func__,
+               mesh_obj->id.name + 2);
+    return;
+  }
+
+  /* Add the armature modifier, if one doesn't exist. */
+  if (!BKE_modifiers_findby_type(mesh_obj, eModifierType_Armature)) {
+    ModifierData *md = BKE_modifier_new(eModifierType_Armature);
+    BLI_addtail(&mesh_obj->modifiers, md);
+  }
+
+  /* Create a deform group per joint. */
+  std::vector<bDeformGroup *> joint_def_grps(joints.size(), nullptr);
+
+  for (int idx : used_indices) {
+    std::string joint_name = pxr::SdfPath(joints[idx]).GetName();
+    if (!BKE_object_defgroup_find_name(mesh_obj, joint_name.c_str())) {
+      bDeformGroup *def_grp = BKE_object_defgroup_add_name(mesh_obj, joint_name.c_str());
+      joint_def_grps[idx] = def_grp;
+    }
+  }
+
+  /* Set the deform group verts and weights. */
+  for (int i = 0; i < mesh->totvert; ++i) {
+    /* Offset into the weights array, which is
+     * always 0 for constant interpolation. */
+    int offset = 0;
+    if (interp == pxr::UsdGeomTokens->vertex) {
+      offset = i * joint_weights_elem_size;
+    }
+    for (int j = 0; j < joint_weights_elem_size; ++j) {
+      const int k = offset + j;
+      const float w = joint_weights[k];
+      if (w < .00001) {
+        /* No deform group if zero weight. */
+        continue;
+      }
+      const int joint_idx = joint_indices[k];
+      if (bDeformGroup *def_grp = joint_def_grps[joint_idx]) {
+        ED_vgroup_vert_add(mesh_obj, def_grp, i, w, WEIGHT_REPLACE);
+      }
+    }
+  }
+}
+
+}  // namespace blender::io::usd
diff --git a/source/blender/io/usd/intern/usd_skel_convert.h b/source/blender/io/usd/intern/usd_skel_convert.h
new file mode 100644
index 00000000000..0cc4b1c3443
--- /dev/null
+++ b/source/blender/io/usd/intern/usd_skel_convert.h
@@ -0,0 +1,72 @@
+/* SPDX-FileCopyrightText: 2023 NVIDIA Corporation. All rights reserved.
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+#pragma once
+
+#include <map>
+#include <pxr/usd/usd/prim.h>
+#include <pxr/usd/usdSkel/skeletonQuery.h>
+
+struct Main;
+struct Object;
+struct Scene;
+struct USDExportParams;
+struct USDImportParams;
+
+namespace blender::io::usd {
+
+struct ImportSettings;
+
+/**
+ * This file contains utilities for converting between UsdSkel data and
+ * Blender armatures and shape keys. The following is a reference on the
+ * UsdSkel API:
+ *
+ * https://openusd.org/23.05/api/usd_skel_page_front.html
+ */
+
+/**
+ * Import USD blend shapes from a USD primitive as shape keys on a mesh
+ * object. Optionally, if the blend shapes have animating weights, the
+ * time-sampled weights will be imported as shape key animation curves.
+ * If the USD primitive does not have blend shape targets defined, this
+ * function is a no-op.
+ *
+ * \param bmain: Main pointer
+ * \param mesh_obj: Mesh object to which imported shape keys will be added
+ * \param prim: The USD primitive from which blendshapes will be imported
+ * \param import_anim: Whether to import time-sampled weights as shape key
+ *                     animation curves
+ */
+void import_blendshapes(Main *bmain,
+                        Object *mesh_obj,
+                        const pxr::UsdPrim &prim,
+                        bool import_anim = true);
+
+/**
+ * Import the given USD skeleton as an armature object. Optionally, if the
+ * skeleton has an animation defined, the time sampled joint transforms will be
+ * imported as bone animation curves.
+ *
+ * \param bmain: Main pointer
+ * \param arm_obj: Armature object to which the bone hierachy will be added
+ * \param skel: The USD skeleton from which bones and animation will be imported
+ * \param import_anim: Whether to import time-sampled joint transforms as bone
+ *                     animation curves
+ */
+void import_skeleton(Main *bmain,
+                     Object *arm_obj,
+                     const pxr::UsdSkelSkeleton &skel,
+                     bool import_anim = true);
+/**
+ * Import skinning data from a source USD prim as deform groups and an armature
+ * modifier on the given mesh object. If the USD prim does not have a skeleton
+ * binding defined, this function is a no-op.
+ *
+ * \param bmain: Main pointer
+ * \param obj: Mesh object to which an armature modifier will be added
+ * \param prim: The USD primitive from which skinning data will be imported
+ */
+void import_mesh_skel_bindings(Main *bmain, Object *mesh_obj, const pxr::UsdPrim &prim);
+
+}  // namespace blender::io::usd
diff --git a/source/blender/io/usd/usd.h b/source/blender/io/usd/usd.h
index 35d32734638..942c6286053 100644
--- a/source/blender/io/usd/usd.h
+++ b/source/blender/io/usd/usd.h
@@ -72,6 +72,8 @@ struct USDImportParams {
   bool import_meshes;
   bool import_volumes;
   bool import_shapes;
+  bool import_skeletons;
+  bool import_blendshapes;
   char *prim_path_mask;
   bool import_subdiv;
   bool import_instance_proxies;
diff --git a/tests/python/bl_usd_import_test.py b/tests/python/bl_usd_import_test.py
index fd183895912..a368fad56fc 100644
--- a/tests/python/bl_usd_import_test.py
+++ b/tests/python/bl_usd_import_test.py
@@ -13,6 +13,8 @@ from pxr import Sdf
 
 import bpy
 
+from mathutils import Matrix, Vector, Quaternion, Euler
+
 args = None
 
 
@@ -302,6 +304,106 @@ class USDImportTest(AbstractUSDTest):
 
         self.assertEqual(4, num_uvmaps_found, "One or more test materials failed to import")
 
+    def test_import_usd_blend_shapes(self):
+        """Test importing USD blend shapes with animated weights."""
+
+        infile = str(self.testdir / "usd_blend_shape_test.usda")
+        res = bpy.ops.wm.usd_import(filepath=infile)
+        self.assertEqual({'FINISHED'}, res)
+
+        obj = bpy.data.objects["Plane"]
+
+        obj.active_shape_key_index = 1
+
+        key = obj.active_shape_key
+        self.assertEqual(key.name, "Key_1", "Unexpected shape key name")
+
+        # Verify the number of shape key points.
+        self.assertEqual(len(key.data), 4, "Unexpected number of shape key point")
+
+        # Verify shape key point coordinates
+
+        # Reference point values.
+        refs = ((-2.51, -1.92, 0.20), (0.86, -1.46, -0.1),
+                (-1.33, 1.29, .84), (1.32, 2.20, -0.42))
+
+        for i in range(4):
+            co = key.data[i].co
+            ref = refs[i]
+            # Compare coordinates.
+            for j in range(3):
+                self.assertAlmostEqual(co[j], ref[j], 2)
+
+        # Verify the shape key values.
+        bpy.context.scene.frame_set(1)
+        self.assertAlmostEqual(key.value, .002, 1)
+        bpy.context.scene.frame_set(30)
+        self.assertAlmostEqual(key.value, .900, 3)
+        bpy.context.scene.frame_set(60)
+        self.assertAlmostEqual(key.value, .100, 3)
+
+    def test_import_usd_skel_joints(self):
+        """Test importing USD animated skeleton joints."""
+
+        infile = str(self.testdir / "arm.usda")
+        res = bpy.ops.wm.usd_import(filepath=infile)
+        self.assertEqual({'FINISHED'}, res)
+
+        # Verify armature was imported.
+        arm_obj = bpy.data.objects["Skel"]
+        self.assertEqual(arm_obj.type, "ARMATURE", "'Skel' object is not an armature")
+
+        arm = arm_obj.data
+        bones = arm.bones
+
+        # Verify bone parenting.
+        self.assertIsNone(bones['Shoulder'].parent, "Shoulder bone should not be parented")
+        self.assertEqual(bones['Shoulder'], bones['Elbow'].parent, "Elbow bone should be child of Shoulder bone")
+        self.assertEqual(bones['Elbow'], bones['Hand'].parent, "Hand bone should be child of Elbow bone")
+
+        # Verify armature modifier was created on the mesh.
+        mesh_obj = bpy.data.objects['Arm']
+        # Get all the armature modifiers on the mesh.
+        arm_mods = [m for m in mesh_obj.modifiers if m.type == "ARMATURE"]
+        self.assertEqual(len(arm_mods), 1, "Didn't get expected armatrue modifier")
+        self.assertEqual(arm_mods[0].object, arm_obj, "Armature modifier does not reference the imported armature")
+
+        # Verify expected deform groups.
+        # There are 4 points in each group.
+        for i in range(4):
+            self.assertAlmostEqual(mesh_obj.vertex_groups['Hand'].weight(i), 1.0, 2, "Unexpected weight for Hand deform vert")
+            self.assertAlmostEqual(mesh_obj.vertex_groups['Shoulder'].weight(4 + i), 1.0, 2, "Unexpected weight for Shoulder deform vert")
+            self.assertAlmostEqual(mesh_obj.vertex_groups['Elbow'].weight(8 + i), 1.0, 2, "Unexpected weight for Elbow deform vert")
+
+        action = bpy.data.actions['SkelAction']
+
+        # Verify the Elbow joint rotation animation.
+        curve_path = 'pose.bones["Elbow"].rotation_quaternion'
+
+        # Quat W
+        f = action.fcurves.find(curve_path, index=0)
+        self.assertIsNotNone(f, "Couldn't find Elbow rotation quaternion W curve")
+        self.assertAlmostEqual(f.evaluate(0), 1.0, 2, "Unexpected value for rotation quaternion W curve at frame 0")
+        self.assertAlmostEqual(f.evaluate(10), 0.707, 2, "Unexpected value for rotation quaternion W curve at frame 10")
+
+        # Quat X
+        f = action.fcurves.find(curve_path, index=1)
+        self.assertIsNotNone(f, "Couldn't find Elbow rotation quaternion X curve")
+        self.assertAlmostEqual(f.evaluate(0), 0.0, 2, "Unexpected value for rotation quaternion X curve at frame 0")
+        self.assertAlmostEqual(f.evaluate(10), 0.707, 2, "Unexpected value for rotation quaternion X curve at frame 10")
+
+        # Quat Y
+        f = action.fcurves.find(curve_path, index=2)
+        self.assertIsNotNone(f, "Couldn't find Elbow rotation quaternion Y curve")
+        self.assertAlmostEqual(f.evaluate(0), 0.0, 2, "Unexpected value for rotation quaternion Y curve at frame 0")
+        self.assertAlmostEqual(f.evaluate(10), 0.0, 2, "Unexpected value for rotation quaternion Y curve at frame 10")
+
+        # Quat Z
+        f = action.fcurves.find(curve_path, index=3)
+        self.assertIsNotNone(f, "Couldn't find Elbow rotation quaternion Z curve")
+        self.assertAlmostEqual(f.evaluate(0), 0.0, 2, "Unexpected value for rotation quaternion Z curve at frame 0")
+        self.assertAlmostEqual(f.evaluate(10), 0.0, 2, "Unexpected value for rotation quaternion Z curve at frame 10")
+
 
 def main():
     global args