An example poseBlending maya plugin using the same idea as my poseblendlib

blendPose.py

import maya.api.OpenMaya as om2
import numpy as np

maya_useNewAPI = True


class blendPose(om2.MPxNode):
    typeName = "blendPose"
    typeId = om2.MTypeId(0x00122717)

    aOutputPose = None
    aOutputRawPose = None
    aOriginalPose = None
    aInputTarget = None
    aInputTargetPose = None
    aInputTargetLevel = None
    aWeight = None

    def getArrayData(self, arrayHandle, getter):
        """Get the data from the given handle as a dictionary of elements"""
        itemCount = len(arrayHandle)
        if itemCount == 0:
            return {}

        poses = {}
        for i in range(itemCount):
            arrayHandle.jumpToPhysicalElement(i)
            key = arrayHandle.elementLogicalIndex()
            dataHandle = arrayHandle.inputValue()
            poses[key] = getter(dataHandle)

        return poses

    def getPoseDict(self, arrayHandle):
        return self.getArrayData(arrayHandle, lambda x: x.asMatrix())

    def getWeights(self, arrayHandle):
        return self.getArrayData(arrayHandle, lambda x: x.asFloat())

    def getAllData(self, plug, dataBlock):
        restArrayHandle = dataBlock.inputArrayValue(self.aOriginalPose)
        restPose = self.getPoseDict(restArrayHandle)

        weightArrayHandle = dataBlock.inputArrayValue(self.aWeight)
        weights = self.getWeights(weightArrayHandle)

        targetArrayHandle = dataBlock.inputArrayValue(self.aInputTarget)
        targetCount = len(targetArrayHandle)
        if targetCount == 0:
            dataBlock.setClean(plug)
            return restPose, {}, {}

        matIdxs = set(restPose.keys())
        targets = {}
        for i in range(targetCount):
            targetArrayHandle.jumpToPhysicalElement(i)
            poseIdx = targetArrayHandle.elementLogicalIndex()
            weight = weights.get(poseIdx, 0.0)
            if abs(weight) <= 1.0e-12:
                continue
            targetIndexHandle = targetArrayHandle.inputValue()
            poseArrayHandle = om2.MArrayDataHandle(
                targetIndexHandle.child(self.aInputTargetPose)
            )

            poseLevel = targetIndexHandle.child(self.aInputTargetLevel).asInt()

            poseMats = self.getPoseDict(poseArrayHandle)
            matIdxs |= poseMats.keys()
            targets[poseIdx] = (poseMats, weight, poseLevel)

        for mi in matIdxs - restPose.keys():
            restPose[mi] = om2.MMatrix()

        return restPose, targets, weights

    def setAllData(
        self,
        plug,
        dataBlock,
        rawMats: dict[int, om2.MMatrix],
        quatMats: dict[int, om2.MMatrix],
    ) -> bool:
        outputArrayHandle = dataBlock.outputArrayValue(self.aOutputPose)
        outputBuilder = outputArrayHandle.builder()

        outputRawArrayHandle = dataBlock.outputArrayValue(self.aOutputRawPose)
        outputRawBuilder = outputRawArrayHandle.builder()

        for idx in rawMats.keys() | quatMats.keys():
            if idx in rawMats:
                outputRawHandle = outputRawBuilder.addElement(idx)
                outputRawHandle.set(rawMats[idx])

            if idx in quatMats:
                outputHandle = outputBuilder.addElement(idx)
                outputHandle.set(quatMats[idx])

        outputRawArrayHandle.setAllClean()
        outputArrayHandle.setAllClean()
        dataBlock.setClean(plug)
        return True

    def computeRawMats(
        self,
        restPose: dict[int, om2.MMatrix],
        targets: dict[int, tuple[dict[int, om2.MMatrix], float, int]],
    ) -> dict[int, om2.MMatrix]:
        ret = {}
        allweight = 0.0
        for pose, weight, _level in targets.values():
            allweight += abs(weight)
            for matIdx, mat in pose.items():
                val = weight * (mat - restPose[matIdx])
                if matIdx in ret:
                    ret[matIdx] += val
                else:
                    ret[matIdx] = val
        if allweight == 0.0:
            return restPose

        for idx in restPose:
            ret[idx] += restPose[idx]

        return ret

    def outerq(self, q):
        out = []
        if q.w < 0:
            q.negateIt()
        for i in q:
            for j in q:
                out.append(i * j)
        return om2.MMatrix(out)

    def computeQuatMats(
        self,
        restPose: dict[int, om2.MMatrix],
        targets: dict[int, tuple[dict[int, om2.MMatrix], float, int]],
    ) -> dict[int, om2.MMatrix]:
        # First build the quaternions in rest-space
        # And sum the weighted outer product (normalized)
        qi = om2.MQuaternion()
        qmats: dict[int, dict[int, om2.MMatrix]] = {}
        levelCounts: dict[int, int] = {}
        for pose, weight, level in targets.values():
            if abs(weight) < 1.0e-12:
                continue
            for matIdx, mat in pose.items():
                matquat = om2.MQuaternion().setValue(mat)
                restquat = om2.MQuaternion().setValue(restPose[matIdx])
                locquat = matquat * restquat.inverse()
                # Gotta do it by the squared weight
                # because we're taking the "outer square"
                outerVal = self.outerq(
                    om2.MQuaternion.slerp(qi, locquat.normal(), weight)
                )

                qmatlevel = qmats.setdefault(level, {})
                if matIdx in qmatlevel:
                    qmatlevel[matIdx] += outerVal
                else:
                    qmatlevel[matIdx] = outerVal
                levelCounts[level] = levelCounts.get(level, 0) + 1

        # Then do the eigen-thing, and slerp to un-normalize the weight
        eigs: dict[int, dict[int, om2.MQuaternion]] = {}
        for level, qmatpose in qmats.items():
            count = levelCounts[level]
            leveleigs = {}
            eigs[level] = leveleigs
            for matIdx, qmat in qmatpose.items():
                wmat = np.array(list(qmat)).reshape((4, 4)).T
                # I don't know why np.linalg.eigh doesn't work
                # The matrices are real and symmetric
                evals, evecs = np.linalg.eig(wmat)
                idx = np.argmax(np.abs(evals))
                q = om2.MQuaternion(*evecs[:, idx])
                q = om2.MQuaternion.slerp(qi, q, count)
                if q.w < 0:
                    q.negateIt()
                leveleigs[matIdx] = q.normal()

        # Then multiply all the layers in order
        qret = {k: om2.MQuaternion().setValue(m) for k, m in restPose.items()}
        for level in sorted(eigs):
            for matIdx, q in eigs[level].items():
                qret[matIdx] = q * qret[matIdx]

        # Finally, convert back to matrices
        return {k: q.asMatrix() for k, q in qret.items()}

    def compute(self, plug, dataBlock):
        alldata = self.getAllData(plug, dataBlock)
        restPose: dict[int, om2.MMatrix] = alldata[0]
        targets: dict[int, tuple[dict[int, om2.MMatrix], float, int]] = alldata[1]
        weights: dict[int, float] = alldata[2]

        if not targets:
            return self.setAllData(plug, dataBlock, restPose, restPose)

        rawMats = self.computeRawMats(restPose, targets)
        quatMats = self.computeQuatMats(restPose, targets)

        for idx in rawMats.keys() & quatMats.keys():
            r = rawMats[idx]
            q = quatMats[idx]
            q[12] = r[12]
            q[13] = r[13]
            q[14] = r[14]

        return self.setAllData(plug, dataBlock, rawMats, quatMats)

    @classmethod
    def creator(cls):
        return blendPose()

    @classmethod
    def initialize(cls):
        nAttr = om2.MFnNumericAttribute()
        mAttr = om2.MFnMatrixAttribute()
        cAttr = om2.MFnCompoundAttribute()

        cls.aOutputPose = mAttr.create(
            "outputPose", "out", om2.MFnMatrixAttribute.kDouble
        )
        mAttr.usesArrayDataBuilder = True
        mAttr.writable = False
        mAttr.array = True
        cls.addAttribute(cls.aOutputPose)

        cls.aOutputRawPose = mAttr.create(
            "outputRawPose", "raw", om2.MFnMatrixAttribute.kDouble
        )
        mAttr.usesArrayDataBuilder = True
        mAttr.writable = False
        mAttr.array = True
        cls.addAttribute(cls.aOutputRawPose)

        cls.aOriginalPose = mAttr.create(
            "originalPose", "og", om2.MFnMatrixAttribute.kDouble
        )
        mAttr.array = True
        mAttr.storable = True
        cls.addAttribute(cls.aOriginalPose)

        cls.aInputTargetPose = mAttr.create(
            "inputTargetPose", "itp", om2.MFnMatrixAttribute.kDouble
        )
        mAttr.array = True
        mAttr.storable = True

        cls.aInputTargetLevel = nAttr.create(
            "inputTargetLevel", "itl", om2.MFnNumericData.kInt, 0
        )
        nAttr.storable = True

        cls.aInputTarget = cAttr.create("inputTarget", "it")
        cAttr.addChild(cls.aInputTargetPose)
        cAttr.addChild(cls.aInputTargetLevel)
        cAttr.array = True
        cls.addAttribute(cls.aInputTarget)

        cls.aWeight = nAttr.create("weight", "w", om2.MFnNumericData.kFloat, 0.0)
        nAttr.setSoftMin(0.0)
        nAttr.setSoftMax(1.0)
        nAttr.array = True
        nAttr.keyable = True
        nAttr.storable = True
        cls.addAttribute(cls.aWeight)

        inputs = [
            cls.aOriginalPose,
            cls.aInputTargetPose,
            cls.aInputTarget,
            cls.aWeight,
        ]
        outputs = [cls.aOutputPose, cls.aOutputRawPose]
        for inat in inputs:
            for outat in outputs:
                cls.attributeAffects(inat, outat)

        return True


def initializePlugin(pluginObj):
    fnPlugin = om2.MFnPlugin(pluginObj, "Tyler Fox", "0.0.1")
    fnPlugin.registerNode(
        blendPose.typeName,
        blendPose.typeId,
        blendPose.creator,
        blendPose.initialize,
        om2.MPxNode.kDependNode,
    )


def uninitializePlugin(pluginObj):
    fnPlugin = om2.MFnPlugin(pluginObj)
    fnPlugin.deregisterNode(blendPose.typeId)