I know that the usual way to perform rotation of an image is to compute the new pixel coordinates by multiplying with a rotation matrix $ \begin{pmatrix}x_{rot} \\ y_{rot} \end{pmatrix} = \begin{pmatrix} cos(a) & -sin(a)\\ sin(a) & cos(a)\end{pmatrix} \begin{pmatrix}x \\ y \end{pmatrix}$ and then, lets say for example for bilinear interpolation, to take the weighted average of the 4 closest pixel intensites (a neighborhood of pixels values located around $\begin{pmatrix}x \\ y \end{pmatrix}$) in the original image (looking "backwards"). Of course here the coordinates x,y, of the pixels are treated distinctly from the intensities at these x,y locations used to perform the weighted averages (i.e., it doesn't make sense to multiply pixel intensities by the rotation matrix nor does it make sense to average pixel coordinates, these are 2 distinct steps). (I put for angles "not a multiple of 90 degrees" since these obviously do not require interpolation).
My question is that I would be interested in a practical implementation (e.g. Matlab or python) of performing interpolation by convolution with an interpolation kernel, much like for image up-sampling (like in here: Why does upsampling and interpolation by convolution introduce a shift compared to imresize?), but for the case of rotations. I've heard/read that this could be done by some kind of matrix multiplication with a matrix filled with the convolution kernel arranged in some (at the moment mysterious way!) for me. I have absolutely no idea and wasn't able to find elsewhere how to do that but I know it can be very useful in some situations. Any ideas? If you have a mathematical formulation going along with the implementation this is very interesting also.
