(Matlab图像旋转(FFT矫正图像倾斜))

% % 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; % %************************************************************************ %* (c) Copyright 1998 * %* Biomathematics Resource * %* Mayo Foundation * %************************************************************************ % % 11/22/98 Implemented by Edward Brian Welch, edwardbrianwelch@yahoo.com % function[Image] = RotateImage(Image, degrees) % NOTE: % This function is partially vectorized (only single loops, % no double loops) to gain some benefit in speed. The tradeoff % is that more memory is required to hold the large matrices. % It would be possible to completely vectorize certain parts % of the mfile, but that would require large amounts of memory % and would also make the code less readable. [xdim ydim] = size(Image); if mod(xdim,2)==1 | mod(ydim,2)==1, error('RotateImage() can only rotate images with even dimensions.'); end % Determine type of image to return if isreal(Image), real_flag = 1; else real_flag = 0; end % Put degrees into the range [0,360] while degrees<0, degrees = degrees + 360; end while degrees>360, degrees = degrees - 360; end % Count number of elementary 90 degree rotations % to put rotation into range [-45,45] count=0; while degrees>45, degrees = degrees - 90; count = count + 1; end Image = rot90(Image, count); % Calculate Trigonometric Values % Not Necessary in Matlab implementation to Negate degrees % so that CCW rotations are positive theta = (degrees/360)*2*pi; sin_theta = sin(theta); negtan_thetadiv2 = -tan(theta/2); %--------------------------- % FIRST SKEW % |x1| = | 1 -tan(a/2) ||x| % |y1| = | 0 1 ||y| %--------------------------- % Pad image rows to double the row size Image2 = zeros(2*xdim,ydim); Image2( (xdim/2+1):(xdim/2+xdim),:)=Image; % Forward FFT image rows Image2 = fft(Image2, 2*xdim, 1); % Calculate image's center coordinates xno = (2*xdim-1)/2; yno = (ydim-1)/2; % Prepare to use the Fourier Shift Theorem % f(x-deltax) <=> exp(-j*2*pi*deltax*k/N)*F(k) % Initialize constant part of the exponent % expression. The skew is row dependent. cons1 = (-2.0*pi/(2*xdim))*negtan_thetadiv2; % Calculate k values (Nyquist is at x=xno) k_array = zeros((2*xdim),1); for x=1:(2*xdim), if (x-1)<=xno, k_array(x) = (x-1); else k_array(x) = (x-1-2*xdim); end end for y=1:ydim, % Skew dependent part of expression cons2 = cons1 * (y-1-yno); % Calculate the angles angle_array = cons2*k_array; % Rotate the complex numbers by those angles sin_ang = sin(angle_array); cos_ang = cos(angle_array); newr = real(Image2(:,y)).*cos_ang - imag(Image2(:,y)).*sin_ang; newi = real(Image2(:,y)).*sin_ang + imag(Image2(:,y)).*cos_ang; Image2(:,y) = newr + newi*i; end %--------------------------- % SECOND SKEW % |x2| = | 1 0 ||x1| % |y2| = | sin(a) 1 ||y1| %--------------------------- % Pad the image columns in hybrid space Image22 = zeros(2*xdim, 2*ydim); Image22(:,(ydim/2+1):(ydim/2+ydim) )=Image2; % Perform a Forward FFT on the image columns Image22 = fft(Image22, 2*ydim, 2); % Mask aliased components from SKEW 1 Image22 = fftshift(Image22); Image22 = MaskImage( Image22, 1, 0, negtan_thetadiv2, 1); Image22 = fftshift(Image22); % Inverse FFT the image rows Image22 = ifft(Image22, 2*xdim, 1); % Calculate image's center coordinates xno = (2*xdim-1)/2; yno = (2*ydim-1)/2; % Initialize constant part of the exponent % expression. The skew is column dependent. cons1 = (-2.0*pi/(2*ydim))*sin_theta; % Calculate k values (Nyquist is at y=yno) k_array = zeros(1,(2*ydim)); for y=1:(2*ydim), if (y-1)<=yno, k_array(y) = (y-1); else k_array(y) = (y-1-2*ydim); end end for x=1:(2*xdim), % Skew dependent part of expression cons2 = cons1 * (x-1-xno); % Calculate the angles angle_array = cons2*k_array; % Rotate the complex numbers by those angles sin_ang = sin(angle_array); cos_ang = cos(angle_array); newr = real(Image22(x,:)).*cos_ang - imag(Image22(x,:)).*sin_ang; newi = real(Image22(x,:)).*sin_ang + imag(Image22(x,:)).*cos_ang; Image22(x,:) = newr + newi*i; end %--------------------------- % THIRD SKEW % |x3| = | 1 -tan(a/2) ||x2| % |y3| = | 0 1 ||y2| %--------------------------- % Forward FFT image rows Image22 = fft(Image22, 2*xdim, 1); % Mask aliased components from SKEW 2 Image22 = fftshift(Image22); Image22 = MaskImage( Image22, 1, sin_theta, 0, 1); Image22 = fftshift(Image22); % Inverse FFT the image columns Image22 = ifft(Image22, 2*ydim, 2); % Crop image columns Image2 = Image22(:,(ydim/2+1):(ydim/2+ydim) ); % Calculate image's center coordinates xno = (2*xdim-1)/2; yno = (ydim-1)/2; % Initialize constant part of the exponent % expression. The skew is row dependent. cons1 = (-2.0*pi/(2*xdim))*negtan_thetadiv2; % Calculate k values (Nyquist is at x=xno) k_array = zeros((2*xdim),1); for x=1:(2*xdim), if (x-1)<=xno, k_array(x) = (x-1); else k_array(x) = (x-1-2*xdim); end end for y=1:ydim, % Skew dependent part of expression cons2 = cons1 * (y-1-yno); % Calculate the angles angle_array = cons2*k_array; % Rotate the complex numbers by those angles sin_ang = sin(angle_array); cos_ang = cos(angle_array); newr = real(Image2(:,y)).*cos_ang - imag(Image2(:,y)).*sin_ang; newi = real(Image2(:,y)).*sin_ang + imag(Image2(:,y)).*cos_ang; Image2(:,y) = newr + newi*i; end % Forward FFT the image columns Image2 = fft(Image2, ydim, 2); % Mask aliased components from SKEW 3 % The /2.0 factor is there because columns have been cropped in half Image2 = fftshift(Image2); Image2 = MaskImage( Image2, 1, 0, negtan_thetadiv2/2.0, 1); Image2 = fftshift(Image2); % Inverse 2D FFT the image Image2 = ifft2(Image2); % Crop the rows to obtain final result Image = Image2( (xdim/2+1):(xdim/2+xdim),:) ; % Return a Real image if original Image was Real if real_flag==1, Image = real(Image); end % % Image_Masked = MaskImage(Image, a, b, c, d) % % Masks an image according to the (x,y) coordinate % transform matrix: |x'| = | a b | * |x| % |y'| | c d | |y| % % (x,y) positions that become (x',y') outside the boundary % of the original image will be set to 0. % % The masking is softened at the edges by a convolution. % % %************************************************************************ %* (c) Copyright 1998 * %*