The correct way doing it is not so trivial.
- The rotation matrix in your post is "centered" - (0,0) is the center coordinate.
We need transformation matrix in which (1,1) is the top left coordinate.
- You need to transform all the coordinates of the image.
In your post you are only transforming coordinates (1,1), (2,2), (3,3)...
- Using forward transformation creates "holes" - not all pixels in the destination image are filled.
We need to use backward transformation, and for every destination pixel, get the coordinate of the source pixel (destination to source transformation is considered "backward transformation").
Using the inverse transformation with transformPointsForward
is equivalent to backward transformation.
Note: Your transformation matrix is actually the "backward transformation".
Here is how to solve it (please read the comments):
I = rgb2gray(imread('peppers.png')); % Read sample image and convert to gray.
[rows, cols] = size(I);
theta = 30;
% The rotation matrix is "centered" - coordinate (0, 0) is the center:
cT = [cosd(theta) -sind(theta) 0; sind(theta) cosd(theta) 0; 0 0 1];
top2cen = [1 0 0
0 1 0
(cols+1)/2 (rows+1)/2 1];
cen2top = [1 0 0
0 1 0
-(cols+1)/2 -(rows+1)/2 1];
% We need the rotation matrix to be "top left" - coordinate (1, 1) is the top left coordinate:
T = cen2top*cT*top2cen; % Note: you don't really need to use matrix multiplications for solving this.
tform = affine2d(T);
% All the combinations of u,v coordinates
[U, V] = meshgrid(1:cols, 1:rows);
% Transform all the (u, v) coordinates of the input I.
[X, Y] = transformPointsForward(tform, U, V);
% Round the coordinates - the interpolation method is going to Nearest Neighbor.
X = round(X);
Y = round(Y);
J = zeros(size(I), 'like', I);
% Limit the X,Y coordinates to the valid range.
limX = max(min(X, cols), 1);
limY = max(min(Y, rows), 1);
% Copy the (u,v) pixel in I to position (x,y) in J.
J(sub2ind(size(I), limY, limX)) = I(sub2ind(size(I), V, U));
% Oops... J has many holes...
figure;imshow(J);
imwrite(J, 'fwJ.png');
% Correct way:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We must inverse the transformation - use backward transformation instead of forward transformation.
%cT = inv([cosd(theta) -sind(theta) 0; sind(theta) cosd(theta) 0; 0 0 1]);
cT = [cosd(theta) sind(theta) 0; -sind(theta) cosd(theta) 0; 0 0 1]; % Inverse transformation matrix.
% Repeate the process with inversed transformation.
T = cen2top*cT*top2cen;
tform = affine2d(T);
[U, V] = meshgrid(1:cols, 1:rows);
% Transform all the (x, y) coordinates of the input I.
[X, Y] = transformPointsForward(tform, U, V); % Name the coordinates U, V
% Round the coordinates - the interpolation method is going to Nearest Neighbor.
X = round(X);
Y = round(Y);
J = zeros(size(I), 'like', I);
% Limit the X,Y coordinates to the valid range.
limX = max(min(X, cols), 1);
limY = max(min(Y, rows), 1);
J(sub2ind(size(I), V, U)) = I(sub2ind(size(I), limY, limX));
% Zero the margins (place zeros where X, Y are outside of the valid range):
J((X < 1) | (Y < 1) | (X > cols) | (Y > rows)) = 0;
figure;imshow(J)
imwrite(J, 'J.png');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Testing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reference:
K = imrotate(I, 30, 'nearest', 'crop');
figure;imshow(K)
% Display the difference from imrotate (images are equal).
figure;imagesc(double(K) - double(J));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Result of forward transformation (the wrong way):

Result of backward transformation (the right way):
