If your transform function is transform(control, image) -> transformed_image, where control is 2 * N the the backward should be:
transform_bwd(grad_transformed_image) -> grad_control, grad_image
If you don’t care about optimizing the image, you may not need to compute grad_image.
The grad_transformed_image will be the same shape as transformed_image and the computed grad_control should be 2 * N (the same shape as control).
One way to compute grad_control would be to first compute the jacobian of the transformation and then multiply by the incoming grad_transformed_image. The jacobian is a matrix of size (2*N) x (M) where M is the number of elements in the transformed image. Each entry in the matrix is a partial derivative: how does that control point affect a specific element in the transformed image. The grad_control would then be the matrix-vector product: jacobian x grad_transformed_image.
In practice, you probably don’t want to explicitly compute the Jacobian, since many of the entires will probably be 0. You can likely simplify the computation if you know properties of your transformation, such as if it’s local (control points only affect nearby pixels).