mriphantom.html

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

0077 nellipses = length(rho);
0078 a   = FOV/2 * E(:,2)';
0079 b   = FOV/2 * E(:,3)';
0080 x0  = FOV/2 * E(:,4)';
0081 y0  = FOV/2 * E(:,5)';
0082 alpha = pi*E(:,6)'/180;
0083 
0084 
0085 <span class="comment">%========================================================================</span>
0086 <span class="comment">% Stretch kspace coordinates to column vector (retain original size info)</span>
0087 <span class="comment">%========================================================================</span>
0088 [mk,nk] = size(kx);
0089 kx = kx(:);
0090 ky = ky(:);
0091 k = kx + j*ky;
0092 klen = length(kx);
0093 
0094 <span class="comment">%========================================================================</span>
0095 <span class="comment">% Calculate angle of the vector to the ellipse center in real space.</span>
0096 <span class="comment">%========================================================================</span>
0097 gamma = atan2(y0, x0);
0098 te = sqrt(x0.^2 + y0.^2);
0099 
0100 <span class="comment">%========================================================================</span>
0101 <span class="comment">% Replicate all the ellipse defining vectors to the number of k space values</span>
0102 <span class="comment">%========================================================================</span>
0103 gamma = gamma(ones(klen,1),:);
0104 te = te(ones(klen,1),:);
0105 alpha = alpha(ones(klen,1),:);
0106 rho = rho(ones(klen,1),:);
0107 a = a(ones(klen,1),:);
0108 b = b(ones(klen,1),:);
0109 
0110 <span class="comment">%========================================================================</span>
0111 <span class="comment">% Calculate k space vectors angle and magnitude</span>
0112 <span class="comment">%========================================================================</span>
0113 kphi = angle(k);
0114 kr = abs(k);
0115 
0116 <span class="comment">%========================================================================</span>
0117 <span class="comment">% Replicate k space and time vectors to the number of ellipses</span>
0118 <span class="comment">%========================================================================</span>
0119 kr = kr(:,ones(1,nellipses));
0120 kphi = kphi(:,ones(1,nellipses));
0121 
0122 <span class="comment">%========================================================================</span>
0123 <span class="comment">% Precalculate some values</span>
0124 <span class="comment">%========================================================================</span>
0125 pi2 = 2*pi;
0126 cose = cos(gamma - kphi);
0127 
0128 <span class="comment">%========================================================================</span>
0129 <span class="comment">% Projection angles:</span>
0130 <span class="comment">% Difference of the k space vector angle and the angles of the ellipses</span>
0131 <span class="comment">%========================================================================</span>
0132 sind = sin(kphi - alpha);
0133 cosd = cos(kphi - alpha);
0134 
0135 akphi = sqrt(a.^2 .* cosd.^2 + b.^2 .* sind.^2);
0136 <span class="comment">%========================================================================</span>
0137 <span class="comment">% Calculate the signals for each ellips and each k space point kx,ky</span>
0138 <span class="comment">% avoid division by zero    besselj(1,0) = 0 so set 0/0 = 1</span>
0139 <span class="comment">%========================================================================</span>
0140 knz = (akphi.*kr ~=0);
0141 bess = ones(size(akphi.*kr));
0142 bess(knz) = besselj(1,pi2.*akphi(knz).*kr(knz)) ./(akphi(knz).*kr(knz));
0143 Sk = exp(-j*pi2.*kr.*te.*cose) .*rho.*a.*b.* bess;
0144 <span class="comment">%========================================================================</span>
0145 <span class="comment">% Sum the signals for all ellipses and reshape the data back to</span>
0146 <span class="comment">% the matrix size of the kx and ky input data</span>
0147 <span class="comment">%========================================================================</span>
0148 Sk = sum(Sk,2);
0149 Sk = reshape(Sk, mk,nk);
0150 <span class="comment">%========================================================================</span>
0151 <span class="comment">% Show the result, including the reconstruction if default Carthesian</span>
0152 <span class="comment">% sampling was used</span>
0153 <span class="comment">%========================================================================</span>
0154 <span class="keyword">if</span>  nargin &gt;= 2
0155     fh = figure;
0156     subplot(1,2,1);
0157     <span class="keyword">if</span> ~isempty(P)
0158         imagesc(abs(P));
0159     <span class="keyword">end</span>;    
0160     title(<span class="string">'Software head phantom'</span>)
0161     axis image
0162     
0163     subplot(1,2,2);
0164     imagesc(abs(Sk));
0165     title(<span class="string">'Simulated raw data'</span>)
0166     <span class="keyword">return</span>
0167 <span class="keyword">else</span>
0168     fh = figure;
0169     subplot(1,3,1);
0170     <span class="keyword">if</span> ~isempty(P)
0171         imagesc(abs(P));
0172     <span class="keyword">end</span>;
0173     axis image
0174 
0175     subplot(1,3,2);
0176     imagesc(abs(Sk));    
0177     title(<span class="string">'Simulated raw data'</span>)
0178     axis image
0179 
0180     subplot(1,3,3);
0181     I = fftshift(fft2(Sk));
0182     imagesc(abs(I));
0183     axis image    
0184     title(<span class="string">'2DFFT Reconstructed image'</span>)
0185 <span class="keyword">end</span>;
0186 <span class="comment">%========================================================================</span>
0187 <span class="comment">% Suppress output if not asked for any</span>
0188 <span class="comment">%========================================================================</span>
0189 <span class="keyword">if</span> nargout &lt; 1
0190     Sk = [];
0191 <span class="keyword">end</span>;
0192 
0193 
0194 <span class="comment">%========================================================================</span>
0195 <span class="comment">% END</span>
0196 <span class="comment">%========================================================================</span>
0197 
0198 <span class="comment">%========================================================================</span>
0199 <span class="comment">% local function Copied from function phantom.m</span>
0200 <span class="comment">% to produce a result without the image toolbox</span>
0201 <span class="comment">%========================================================================</span>
0202 
0203 <span class="comment">%========================================================================</span>
0204 <a name="_sub1" href="#_subfunctions" class="code">function toft=modified_shepp_logan</a>
0205 <span class="comment">%</span>
0206 <span class="comment">%   This head phantom is the same as the Shepp-Logan except</span>
0207 <span class="comment">%   the intensities are changed to yield higher contrast in</span>
0208 <span class="comment">%   the image.  Taken from Toft, 199-200.</span>
0209 <span class="comment">%</span>
0210 <span class="comment">%         A    a     b    x0    y0    phi</span>
0211 <span class="comment">%        ---------------------------------</span>
0212 toft = [  1   .69   .92    0     0     0   
0213         -.8  .6624 .8740   0  -.0184   0
0214         -.2  .1100 .3100  .22    0    -18
0215         -.2  .1600 .4100 -.22    0     18
0216          .1  .2100 .2500   0    .35    0
0217          .1  .0460 .0460   0    .1     0
0218          .1  .0460 .0460   0   -.1     0
0219          .1  .0460 .0230 -.08  -.605   0 
0220          .1  .0230 .0230   0   -.606   0
0221          .1  .0230 .0460  .06  -.605   0   ];
0222 <span class="comment">%========================================================================</span>
0223 <span class="comment">% end local function toft=modified_shepp_logan</span>
0224 <span class="comment">%========================================================================</span></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>