@Tripp, I have found that bug too, I think it suffice to replace M=[] with M=eye(4).
There is however another bug I am afraid. I have not gone deep into the code yet but it seems that, when doing 3D NonRigid registration only 1 direction out of 3 is computed, i.e. the resulting grid only has the estimated deformation in the third component, and components 1 and 2 remain unchanged. Ani hint on how to fix this?