RotateImage - Rotates an image by X degrees Image_Rotated = RotateImage(Image, degrees) Rotates an image by the angle degrees in the CCW direction. Degrees may be any number. The function will put degrees in the range 0 to 360 degrees and then into a range of -45 to 45 degrees after performing elementary 90 degree rotations. The rotation performed is the 3 Pass Separable rotation using FFT-based methods to perform the skews. Aliased spectral components are masked after EACH skew using the MaskImage function. Normal Rotation |x'| = | cos(a) -sin(a) ||x| |y'| | sin(a) cos(a) ||y| 3 Pass Rotation |x'| = | 1 -tan(a/2) || 1 0 || 1 -tan(a/2) ||x| |y'| | 0 1 || sin(a) 1 || 0 1 ||y| This function pads images during the rotation process. Images can be non-square, but odd dimensions will cause incorrect padding Therefore an error is returned, if a dimension is odd. Dimensions need not be powers of 2 and will not necessarily be padded to a power of 2. The Matlab fft functions are capable of non-power of 2 transforms. Note: The results often display some Gibbs ringing. Example: load mri img = double( squeeze(D(:,:,1,16)) ); img_rot = real( RotateImage(img,45) ); figure; subplot(1,2,1); imagesc(img,[0 88]); axis image; subplot(1,2,2); imagesc(img_rot,[0 88]); axis image;
0001 function[Image] = RotateImage(Image, degrees) 0002 % RotateImage - Rotates an image by X degrees 0003 % 0004 % Image_Rotated = RotateImage(Image, degrees) 0005 % 0006 % Rotates an image by the angle degrees in the 0007 % CCW direction. Degrees may be any number. 0008 % The function will put degrees in the range 0 0009 % to 360 degrees and then into a range of -45 to 45 0010 % degrees after performing elementary 90 degree rotations. 0011 % 0012 % The rotation performed is the 3 Pass Separable rotation 0013 % using FFT-based methods to perform the skews. Aliased spectral 0014 % components are masked after EACH skew using the MaskImage function. 0015 % 0016 % Normal Rotation 0017 % 0018 % |x'| = | cos(a) -sin(a) ||x| 0019 % |y'| | sin(a) cos(a) ||y| 0020 % 0021 % 3 Pass Rotation 0022 % 0023 % |x'| = | 1 -tan(a/2) || 1 0 || 1 -tan(a/2) ||x| 0024 % |y'| | 0 1 || sin(a) 1 || 0 1 ||y| 0025 % 0026 % This function pads images during the rotation process. Images 0027 % can be non-square, but odd dimensions will cause incorrect padding 0028 % Therefore an error is returned, if a dimension is odd. 0029 % 0030 % Dimensions need not be powers of 2 and will not necessarily be padded 0031 % to a power of 2. The Matlab fft functions are capable of non-power of 2 0032 % transforms. 0033 % 0034 % Note: The results often display some Gibbs ringing. 0035 % 0036 % Example: 0037 % 0038 % load mri 0039 % img = double( squeeze(D(:,:,1,16)) ); 0040 % img_rot = real( RotateImage(img,45) ); 0041 % figure; 0042 % subplot(1,2,1); imagesc(img,[0 88]); axis image; 0043 % subplot(1,2,2); imagesc(img_rot,[0 88]); axis image; 0044 0045 % 0046 %************************************************************************ 0047 %* (c) Copyright 1998 * 0048 %* Biomathematics Resource * 0049 %* Mayo Foundation * 0050 %************************************************************************ 0051 % 0052 % 11/22/98 Implemented by Edward Brian Welch, edwardbrianwelch@yahoo.com 0053 % 0054 0055 0056 % NOTE: 0057 % This function is partially vectorized (only single loops, 0058 % no double loops) to gain some benefit in speed. The tradeoff 0059 % is that more memory is required to hold the large matrices. 0060 % It would be possible to completely vectorize certain parts 0061 % of the mfile, but that would require large amounts of memory 0062 % and would also make the code less readable. 0063 0064 [xdim ydim] = size(Image); 0065 0066 if mod(xdim,2)==1 | mod(ydim,2)==1, 0067 error('RotateImage() can only rotate images with even dimensions.'); 0068 end 0069 0070 % Determine type of image to return 0071 if isreal(Image), 0072 real_flag = 1; 0073 else 0074 real_flag = 0; 0075 end 0076 0077 % Put degrees into the range [0,360] 0078 while degrees<0, 0079 degrees = degrees + 360; 0080 end 0081 0082 while degrees>360, 0083 degrees = degrees - 360; 0084 end 0085 0086 % Count number of elementary 90 degree rotations 0087 % to put rotation into range [-45,45] 0088 count=0; 0089 while degrees>45, 0090 degrees = degrees - 90; 0091 count = count + 1; 0092 end 0093 0094 Image = rot90(Image, count); 0095 0096 % Calculate Trigonometric Values 0097 % Not Necessary in Matlab implementation to Negate degrees 0098 % so that CCW rotations are positive 0099 theta = (degrees/360)*2*pi; 0100 sin_theta = sin(theta); 0101 negtan_thetadiv2 = -tan(theta/2); 0102 0103 %--------------------------- 0104 % FIRST SKEW 0105 % |x1| = | 1 -tan(a/2) ||x| 0106 % |y1| = | 0 1 ||y| 0107 %--------------------------- 0108 0109 % Pad image rows to double the row size 0110 Image2 = zeros(2*xdim,ydim); 0111 Image2( (xdim/2+1):(xdim/2+xdim),:)=Image; 0112 0113 % Forward FFT image rows 0114 Image2 = fft(Image2, 2*xdim, 1); 0115 0116 % Calculate image's center coordinates 0117 xno = (2*xdim-1)/2; 0118 yno = (ydim-1)/2; 0119 0120 % Prepare to use the Fourier Shift Theorem 0121 % f(x-deltax) <=> exp(-j*2*pi*deltax*k/N)*F(k) 0122 0123 % Initialize constant part of the exponent 0124 % expression. The skew is row dependent. 0125 cons1 = (-2.0*pi/(2*xdim))*negtan_thetadiv2; 0126 0127 % Calculate k values (Nyquist is at x=xno) 0128 k_array = zeros((2*xdim),1); 0129 0130 for x=1:(2*xdim), 0131 if (x-1)<=xno, 0132 k_array(x) = (x-1); 0133 else 0134 k_array(x) = (x-1-2*xdim); 0135 end 0136 end 0137 0138 for y=1:ydim, 0139 % Skew dependent part of expression 0140 cons2 = cons1 * (y-1-yno); 0141 0142 % Calculate the angles 0143 angle_array = cons2*k_array; 0144 0145 % Rotate the complex numbers by those angles 0146 sin_ang = sin(angle_array); 0147 cos_ang = cos(angle_array); 0148 newr = real(Image2(:,y)).*cos_ang - imag(Image2(:,y)).*sin_ang; 0149 newi = real(Image2(:,y)).*sin_ang + imag(Image2(:,y)).*cos_ang; 0150 Image2(:,y) = newr + newi*i; 0151 end 0152 0153 %--------------------------- 0154 % SECOND SKEW 0155 % |x2| = | 1 0 ||x1| 0156 % |y2| = | sin(a) 1 ||y1| 0157 %--------------------------- 0158 0159 % Pad the image columns in hybrid space 0160 Image22 = zeros(2*xdim, 2*ydim); 0161 Image22(:,(ydim/2+1):(ydim/2+ydim) )=Image2; 0162 0163 % Perform a Forward FFT on the image columns 0164 Image22 = fft(Image22, 2*ydim, 2); 0165 0166 % Mask aliased components from SKEW 1 0167 Image22 = fftshift(Image22); 0168 Image22 = MaskImage( Image22, 1, 0, negtan_thetadiv2, 1); 0169 Image22 = fftshift(Image22); 0170 0171 % Inverse FFT the image rows 0172 Image22 = ifft(Image22, 2*xdim, 1); 0173 0174 % Calculate image's center coordinates 0175 xno = (2*xdim-1)/2; 0176 yno = (2*ydim-1)/2; 0177 0178 % Initialize constant part of the exponent 0179 % expression. The skew is column dependent. 0180 cons1 = (-2.0*pi/(2*ydim))*sin_theta; 0181 0182 % Calculate k values (Nyquist is at y=yno) 0183 k_array = zeros(1,(2*ydim)); 0184 for y=1:(2*ydim), 0185 if (y-1)<=yno, 0186 k_array(y) = (y-1); 0187 else 0188 k_array(y) = (y-1-2*ydim); 0189 end 0190 end 0191 0192 for x=1:(2*xdim), 0193 % Skew dependent part of expression 0194 cons2 = cons1 * (x-1-xno); 0195 0196 % Calculate the angles 0197 angle_array = cons2*k_array; 0198 0199 % Rotate the complex numbers by those angles 0200 sin_ang = sin(angle_array); 0201 cos_ang = cos(angle_array); 0202 newr = real(Image22(x,:)).*cos_ang - imag(Image22(x,:)).*sin_ang; 0203 newi = real(Image22(x,:)).*sin_ang + imag(Image22(x,:)).*cos_ang; 0204 Image22(x,:) = newr + newi*i; 0205 end 0206 0207 %--------------------------- 0208 % THIRD SKEW 0209 % |x3| = | 1 -tan(a/2) ||x2| 0210 % |y3| = | 0 1 ||y2| 0211 %--------------------------- 0212 0213 % Forward FFT image rows 0214 Image22 = fft(Image22, 2*xdim, 1); 0215 0216 % Mask aliased components from SKEW 2 0217 Image22 = fftshift(Image22); 0218 Image22 = MaskImage( Image22, 1, sin_theta, 0, 1); 0219 Image22 = fftshift(Image22); 0220 0221 % Inverse FFT the image columns 0222 Image22 = ifft(Image22, 2*ydim, 2); 0223 0224 % Crop image columns 0225 Image2 = Image22(:,(ydim/2+1):(ydim/2+ydim) ); 0226 0227 % Calculate image's center coordinates 0228 xno = (2*xdim-1)/2; 0229 yno = (ydim-1)/2; 0230 0231 % Initialize constant part of the exponent 0232 % expression. The skew is row dependent. 0233 cons1 = (-2.0*pi/(2*xdim))*negtan_thetadiv2; 0234 0235 % Calculate k values (Nyquist is at x=xno) 0236 k_array = zeros((2*xdim),1); 0237 0238 for x=1:(2*xdim), 0239 0240 if (x-1)<=xno, 0241 k_array(x) = (x-1); 0242 else 0243 k_array(x) = (x-1-2*xdim); 0244 end 0245 end 0246 0247 for y=1:ydim, 0248 % Skew dependent part of expression 0249 cons2 = cons1 * (y-1-yno); 0250 0251 % Calculate the angles 0252 angle_array = cons2*k_array; 0253 0254 % Rotate the complex numbers by those angles 0255 sin_ang = sin(angle_array); 0256 cos_ang = cos(angle_array); 0257 newr = real(Image2(:,y)).*cos_ang - imag(Image2(:,y)).*sin_ang; 0258 newi = real(Image2(:,y)).*sin_ang + imag(Image2(:,y)).*cos_ang; 0259 Image2(:,y) = newr + newi*i; 0260 end 0261 0262 % Forward FFT the image columns 0263 Image2 = fft(Image2, ydim, 2); 0264 0265 % Mask aliased components from SKEW 3 0266 % The /2.0 factor is there because columns have been cropped in half 0267 Image2 = fftshift(Image2); 0268 Image2 = MaskImage( Image2, 1, 0, negtan_thetadiv2/2.0, 1); 0269 Image2 = fftshift(Image2); 0270 0271 % Inverse 2D FFT the image 0272 Image2 = ifft2(Image2); 0273 0274 % Crop the rows to obtain final result 0275 Image = Image2( (xdim/2+1):(xdim/2+xdim),:) ; 0276 0277 % Return a Real image if original Image was Real 0278 if real_flag==1, 0279 Image = real(Image); 0280 end 0281 0282 % 0283 % Image_Masked = MaskImage(Image, a, b, c, d) 0284 % 0285 % Masks an image according to the (x,y) coordinate 0286 % transform matrix: |x'| = | a b | * |x| 0287 % |y'| | c d | |y| 0288 % 0289 % (x,y) positions that become (x',y') outside the boundary 0290 % of the original image will be set to 0. 0291 % 0292 % The masking is softened at the edges by a convolution. 0293 % 0294 0295 % 0296 %************************************************************************ 0297 %* (c) Copyright 1998 * 0298 %* Biomathematics Resource * 0299 %* Mayo Foundation * 0300 %************************************************************************ 0301 % 0302 % 11/22/98 Implemented by Edward Brian Welch, edwardbrianwelch@yahoo.com 0303 % 0304 function[Image] = MaskImage(Image, a, b, c, d) 0305 0306 % NOTE: 0307 % This function is highly vectorized (no loops) 0308 % to gain some benefit in speed. The tradeoff is that 0309 % more memory is required to hold the large matrices 0310 0311 [xdim ydim] = size(Image); 0312 0313 % Calculate center of the image 0314 xno = (xdim-1)/2; 0315 yno = (ydim-1)/2; 0316 0317 x_values = 1:1:xdim; 0318 x_values = x_values - 1 - xno; 0319 x_mat = repmat(x_values',1,ydim); 0320 0321 y_values = 1:1:ydim; 0322 y_values = y_values - 1 - yno; 0323 y_mat = repmat(y_values,xdim,1); 0324 0325 % Calculate new x values for this column 0326 new_x_mat = a*x_mat + b*y_mat; 0327 0328 % Calculate new y values for this column 0329 new_y_mat = c*x_mat + d*y_mat; 0330 0331 % Points whos x or y position is greater than 0332 % absolute value of xno or yno respectively 0333 % should be masked (set to zero). 0334 new_x_mat = abs(new_x_mat) - xno; 0335 new_y_mat = abs(new_y_mat) - yno; 0336 0337 % Only valid points will have a non-zero value after 0338 % the operation below 0339 new_x_mat = sign(new_x_mat)-1; 0340 new_y_mat = sign(new_y_mat)-1; 0341 0342 % ANDing the matrices shows good points 0343 mask = new_x_mat & new_y_mat; 0344 0345 % Soften mask edges 0346 Nx=ceil(xdim*.2); 0347 Ny=ceil(ydim*.2); 0348 mask = conv2(mask,ones(Nx,Ny)/(Nx*Ny),'same'); 0349 0350 Image = Image.*mask; 0351 0352 0353