Nice work
Can we still use WreathProduct's here?
Yes. The permutations act on 4*N points instead of N^4.
WreathProduct(SymmetricGroup(N), P) works ... with a few options for the permutation group, P.
P = Group((1,3)(2,4), (1,4))
P = Group((1,4)(2,3), (1,2))
For irregular box sizes ... (Nr x Nc), you can use
DirectProduct(WreathProduct(SymmetricGroup(Nr), SymmetricGroup(2)), WreathProduct(SymmetricGroup(Nc), SymmetricGroup(2)));
With Nr=Nc=N, that would correspond to P = Group((1,2), (3,4)) ... with ReflectMR and ReflectMC, but no Rotate and no Flip.