Calculating Bending Energy regularization in 2D

Hi, I'm doing a non-rigid registration using B-Splines and I'd like to add a bending energy regularization term in my cost function. Basically this regularization term consists in implementing this equation: , where T is the matrix of displacements and V is the size of the volume. I've done a implementation in 2D, but it not seems correct to me. So how to implement this correctly in 2D. Thanks for the support.

    %Calculating the penalty terms
    dx=[1;0;-1];
    dy=dx.';
    dIdxdy=convn(  convn(grid,dx,'same')   ,dy,'same');
    dIdx2=convn(convn(grid, dx, 'same'), dx, 'same');
    dIdy2=convn(convn(grid, dx, 'same'), dx, 'same');
    
    C = (1/I)*sum(sum(sum(dIdx2.^2+dIdy2.^2+2*dIdxdy.^2)));