rotateimage.html

mri_toolbox是一个工具用来MRI. 来自于SourceForge, 我上传这个软件,希望能结识对医疗软件感兴趣的兄弟.

0136 <span class="keyword">end</span>
0137 
0138 <span class="keyword">for</span> y=1:ydim,
0139    <span class="comment">% Skew dependent part of expression</span>
0140     cons2 = cons1 * (y-1-yno);
0141    
0142    <span class="comment">% Calculate the angles</span>
0143    angle_array = cons2*k_array;
0144    
0145    <span class="comment">% Rotate the complex numbers by those angles</span>
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 <span class="keyword">end</span>
0152 
0153 <span class="comment">%---------------------------</span>
0154 <span class="comment">% SECOND SKEW</span>
0155 <span class="comment">% |x2| = | 1       0 ||x1|</span>
0156 <span class="comment">% |y2| = | sin(a)  1 ||y1|</span>
0157 <span class="comment">%---------------------------</span>
0158 
0159 <span class="comment">% Pad the image columns in hybrid space</span>
0160 Image22 = zeros(2*xdim, 2*ydim);
0161 Image22(:,(ydim/2+1):(ydim/2+ydim) )=Image2;
0162 
0163 <span class="comment">% Perform a Forward FFT on the image columns</span>
0164 Image22 = fft(Image22, 2*ydim, 2);
0165 
0166 <span class="comment">% Mask aliased components from SKEW 1</span>
0167 Image22 = fftshift(Image22);
0168 Image22 = <a href="#_sub1" class="code" title="subfunction[Image] = MaskImage(Image, a, b, c, d)">MaskImage</a>( Image22, 1, 0, negtan_thetadiv2, 1);
0169 Image22 = fftshift(Image22);
0170 
0171 <span class="comment">%  Inverse FFT the image rows</span>
0172 Image22 = ifft(Image22, 2*xdim, 1);
0173 
0174 <span class="comment">% Calculate image's center coordinates</span>
0175 xno = (2*xdim-1)/2;
0176 yno = (2*ydim-1)/2;
0177 
0178 <span class="comment">% Initialize constant part of the exponent</span>
0179 <span class="comment">% expression.  The skew is column dependent.</span>
0180 cons1 = (-2.0*pi/(2*ydim))*sin_theta;
0181 
0182 <span class="comment">% Calculate k values (Nyquist is at y=yno)</span>
0183 k_array = zeros(1,(2*ydim));
0184 <span class="keyword">for</span> y=1:(2*ydim),
0185       <span class="keyword">if</span> (y-1)&lt;=yno,
0186       k_array(y) = (y-1);
0187    <span class="keyword">else</span>
0188       k_array(y) = (y-1-2*ydim);
0189    <span class="keyword">end</span>
0190 <span class="keyword">end</span>
0191 
0192 <span class="keyword">for</span> x=1:(2*xdim),  
0193    <span class="comment">% Skew dependent part of expression</span>
0194    cons2 = cons1 * (x-1-xno);
0195    
0196    <span class="comment">% Calculate the angles</span>
0197    angle_array = cons2*k_array;
0198    
0199    <span class="comment">% Rotate the complex numbers by those angles</span>
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 <span class="keyword">end</span>
0206 
0207 <span class="comment">%---------------------------</span>
0208 <span class="comment">% THIRD SKEW</span>
0209 <span class="comment">% |x3| = | 1  -tan(a/2) ||x2|</span>
0210 <span class="comment">% |y3| = | 0     1      ||y2|</span>
0211 <span class="comment">%---------------------------</span>
0212 
0213 <span class="comment">% Forward FFT image rows</span>
0214 Image22 = fft(Image22, 2*xdim, 1);
0215 
0216 <span class="comment">% Mask aliased components from SKEW 2</span>
0217 Image22 = fftshift(Image22);
0218 Image22 = <a href="#_sub1" class="code" title="subfunction[Image] = MaskImage(Image, a, b, c, d)">MaskImage</a>( Image22, 1, sin_theta, 0,  1);
0219 Image22 = fftshift(Image22);
0220 
0221 <span class="comment">% Inverse FFT the image columns</span>
0222 Image22 = ifft(Image22, 2*ydim, 2);
0223 
0224 <span class="comment">% Crop image columns</span>
0225 Image2 = Image22(:,(ydim/2+1):(ydim/2+ydim) );
0226 
0227 <span class="comment">% Calculate image's center coordinates</span>
0228 xno = (2*xdim-1)/2;
0229 yno = (ydim-1)/2;
0230 
0231 <span class="comment">% Initialize constant part of the exponent</span>
0232 <span class="comment">% expression.  The skew is row dependent.</span>
0233 cons1 = (-2.0*pi/(2*xdim))*negtan_thetadiv2;
0234 
0235 <span class="comment">% Calculate k values (Nyquist is at x=xno)</span>
0236 k_array = zeros((2*xdim),1);
0237 
0238 <span class="keyword">for</span> x=1:(2*xdim),
0239       
0240       <span class="keyword">if</span> (x-1)&lt;=xno,
0241          k_array(x) = (x-1);
0242       <span class="keyword">else</span>
0243          k_array(x) = (x-1-2*xdim);
0244       <span class="keyword">end</span>
0245 <span class="keyword">end</span>
0246 
0247 <span class="keyword">for</span> y=1:ydim,
0248    <span class="comment">% Skew dependent part of expression</span>
0249     cons2 = cons1 * (y-1-yno);
0250    
0251    <span class="comment">% Calculate the angles</span>
0252    angle_array = cons2*k_array;
0253    
0254    <span class="comment">% Rotate the complex numbers by those angles</span>
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 <span class="keyword">end</span>
0261 
0262 <span class="comment">% Forward FFT the image columns</span>
0263 Image2 = fft(Image2, ydim, 2);
0264 
0265 <span class="comment">% Mask aliased components from SKEW 3</span>
0266 <span class="comment">% The /2.0 factor is there because columns have been cropped in half</span>
0267 Image2 = fftshift(Image2);
0268 Image2 = <a href="#_sub1" class="code" title="subfunction[Image] = MaskImage(Image, a, b, c, d)">MaskImage</a>( Image2, 1, 0, negtan_thetadiv2/2.0, 1);
0269 Image2 = fftshift(Image2);
0270 
0271 <span class="comment">% Inverse 2D FFT the image</span>
0272 Image2 = ifft2(Image2);
0273 
0274 <span class="comment">% Crop the rows to obtain final result</span>
0275 Image = Image2( (xdim/2+1):(xdim/2+xdim),:) ;
0276 
0277 <span class="comment">% Return a Real image if original Image was Real</span>
0278 <span class="keyword">if</span> real_flag==1,
0279    Image = real(Image);
0280 <span class="keyword">end</span>
0281 
0282 <span class="comment">%</span>
0283 <span class="comment">%  Image_Masked = MaskImage(Image, a, b, c, d)</span>
0284 <span class="comment">%</span>
0285 <span class="comment">%  Masks an image according to the (x,y) coordinate</span>
0286 <span class="comment">%  transform matrix:  |x'| = | a  b | * |x|</span>
0287 <span class="comment">%                     |y'|   | c  d |   |y|</span>
0288 <span class="comment">%</span>
0289 <span class="comment">%  (x,y) positions that become (x',y') outside the boundary</span>
0290 <span class="comment">%  of the original image will be set to 0.</span>
0291 <span class="comment">%</span>
0292 <span class="comment">%  The masking is softened at the edges by a convolution.</span>
0293 <span class="comment">%</span>
0294 
0295 <span class="comment">%</span>
0296 <span class="comment">%************************************************************************</span>
0297 <span class="comment">%*                    (c) Copyright 1998                                *</span>
0298 <span class="comment">%*                  Biomathematics Resource                                *</span>
0299 <span class="comment">%*                      Mayo Foundation                                  *</span>
0300 <span class="comment">%************************************************************************</span>
0301 <span class="comment">%</span>
0302 <span class="comment">% 11/22/98  Implemented by Edward Brian Welch, [e-mail]</span>
0303 <span class="comment">%</span>
0304 <a name="_sub1" href="#_subfunctions" class="code">function[Image] = MaskImage(Image, a, b, c, d)</a>
0305 
0306 <span class="comment">% NOTE:</span>
0307 <span class="comment">% This function is highly vectorized (no loops)</span>
0308 <span class="comment">% to gain some benefit in speed.  The tradeoff is that</span>
0309 <span class="comment">% more memory is required to hold the large matrices</span>
0310 
0311 [xdim ydim] = size(Image);
0312 
0313 <span class="comment">% Calculate center of the image</span>
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 <span class="comment">% Calculate new x values for this column</span>
0326 new_x_mat = a*x_mat + b*y_mat;
0327 
0328 <span class="comment">% Calculate new y values for this column</span>
0329 new_y_mat = c*x_mat + d*y_mat;
0330 
0331 <span class="comment">% Points whos x or y position is greater than</span>
0332 <span class="comment">% absolute value of xno or yno respectively</span>
0333 <span class="comment">% should be masked (set to zero).</span>
0334 new_x_mat = abs(new_x_mat) - xno;
0335 new_y_mat = abs(new_y_mat) - yno;
0336 
0337 <span class="comment">% Only valid points will have a non-zero value after</span>
0338 <span class="comment">% the operation below</span>
0339 new_x_mat = sign(new_x_mat)-1;
0340 new_y_mat = sign(new_y_mat)-1;
0341 
0342 <span class="comment">% ANDing the matrices shows good points</span>
0343 mask = new_x_mat &amp; new_y_mat;
0344 
0345 <span class="comment">% Soften mask edges</span>
0346 Nx=ceil(xdim*.2);
0347 Ny=ceil(ydim*.2);
0348 mask = conv2(mask,ones(Nx,Ny)/(Nx*Ny),<span class="string">'same'</span>);
0349 
0350 Image = Image.*mask;
0351   
0352 
0353</pre></div>
<hr><address>Generated on Fri 21-May-2004 12:38:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>