MRIPHANTOM - simulated MRI of a Shepp Logan head phantom
Sk = mriphantom(kx,ky)
Creates simulated raw MRI data of a Shepp Logan head phantom
for given k space points
input is the kx and ky coordinates in the image frequency space
default values for kx and ky are for a Carthesian sampling.
Input
kx, ky Matrices with kx and ky sampling coordinates
E Optional analytical phantom data matrix as produced by
image toolbox function phantom.m (see comments below).
Data are from an analytical expression for a continuous phantom.
Core equations for this are taken form :
Rik van der Walle et al IEEE Trans Med Imag 19(12) p1160.
Note that the standard matlab function just gives the phantom
sampled in image space on a carthesian grid.
This function delivers the simulated raw MRI data sampled through
a user defined path in k-space (image frequency space)
Author : Ronald Ouwekerk Johns Hopkins University Dept. Radiology
601 N. Caroline Street Room 4250 Baltimore, MD 21287-0845 USA
rouwerke@mri.jhu.edu
Written :May 2002
Modified June 2002 : Revised to give a non-NaN result for k=0
Suppress output if nargout == 0
Added comments to explain the difference with phantom.m
Whenever this software is used for publications an acknowledgment
would be much appreciated.
RO0001 function Sk = mriphantom(kx,ky) 0002 0003 % MRIPHANTOM - simulated MRI of a Shepp Logan head phantom 0004 % 0005 % Sk = mriphantom(kx,ky) 0006 % 0007 % Creates simulated raw MRI data of a Shepp Logan head phantom 0008 % for given k space points 0009 % input is the kx and ky coordinates in the image frequency space 0010 % default values for kx and ky are for a Carthesian sampling. 0011 % Input 0012 % kx, ky Matrices with kx and ky sampling coordinates 0013 % E Optional analytical phantom data matrix as produced by 0014 % image toolbox function phantom.m (see comments below). 0015 % Data are from an analytical expression for a continuous phantom. 0016 % Core equations for this are taken form : 0017 % Rik van der Walle et al IEEE Trans Med Imag 19(12) p1160. 0018 % 0019 % Note that the standard matlab function just gives the phantom 0020 % sampled in image space on a carthesian grid. 0021 % This function delivers the simulated raw MRI data sampled through 0022 % a user defined path in k-space (image frequency space) 0023 % 0024 % Author : Ronald Ouwekerk Johns Hopkins University Dept. Radiology 0025 % 601 N. Caroline Street Room 4250 Baltimore, MD 21287-0845 USA 0026 % rouwerke@mri.jhu.edu 0027 % Written :May 2002 0028 % Modified June 2002 : Revised to give a non-NaN result for k=0 0029 % Suppress output if nargout == 0 0030 % Added comments to explain the difference with phantom.m 0031 % 0032 % Whenever this software is used for publications an acknowledgment 0033 % would be much appreciated. 0034 % RO 0035 0036 FOV = 2; 0037 N = 64; 0038 0039 if nargin < 2 0040 %===================================================================== 0041 % Set the field of view to 24 cm (see image toolbox phantom.m) 0042 % resolution 128x128 points 0043 % create k-space coordinates for Carthesion sampling 0044 % Should be equivalent to IFFT of the phantom 0045 %===================================================================== 0046 FOV = 0.24; 0047 res = FOV/N; 0048 Kmax = 1/(2*res); 0049 ky = linspace(-Kmax,Kmax, N); 0050 ky = ky(ones(N,1),:); 0051 kx = permute(ky, [2,1]); 0052 end; 0053 if ( (nargin < 3 ) | isempty(E) ) & (exist('phantom.m')==2) 0054 %======================================================================== 0055 % Create 2D phantom to get defining matrix E 0056 %======================================================================== 0057 [P,E] = phantom; 0058 elseif ~(exist('phantom.m')==2) 0059 P = []; 0060 E = modified_shepp_logan; 0061 end; 0062 0063 %======================================================================== 0064 % copy the elements of E and scale to field of view 0065 % Matlab phantom.m default yields FOV =[-1,1]= 2 0066 % Scale all dimensions to desired FOV (*FOV/2). 0067 % Column 1: rho the additive intensity value of the ellipse 0068 % Column 2: a the length of the horizontal semi-axis of the ellipse 0069 % Column 3: b the length of the vertical semi-axis of the ellipse 0070 % Column 4: x0 the x-coordinate of the center of the ellipse 0071 % Column 5: y0 the y-coordinate of the center of the ellipse 0072 % Column 6: alpha the angle (in degrees) between the horizontal semi-axis 0073 % of the ellipse and the x-axis of the image 0074 %======================================================================== 0075 [me,ne] = size(E); 0076 rho = E(:,1)'; 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 %======================================================================== 0086 % Stretch kspace coordinates to column vector (retain original size info) 0087 %======================================================================== 0088 [mk,nk] = size(kx); 0089 kx = kx(:); 0090 ky = ky(:); 0091 k = kx + j*ky; 0092 klen = length(kx); 0093 0094 %======================================================================== 0095 % Calculate angle of the vector to the ellipse center in real space. 0096 %======================================================================== 0097 gamma = atan2(y0, x0); 0098 te = sqrt(x0.^2 + y0.^2); 0099 0100 %======================================================================== 0101 % Replicate all the ellipse defining vectors to the number of k space values 0102 %======================================================================== 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 %======================================================================== 0111 % Calculate k space vectors angle and magnitude 0112 %======================================================================== 0113 kphi = angle(k); 0114 kr = abs(k); 0115 0116 %======================================================================== 0117 % Replicate k space and time vectors to the number of ellipses 0118 %======================================================================== 0119 kr = kr(:,ones(1,nellipses)); 0120 kphi = kphi(:,ones(1,nellipses)); 0121 0122 %======================================================================== 0123 % Precalculate some values 0124 %======================================================================== 0125 pi2 = 2*pi; 0126 cose = cos(gamma - kphi); 0127 0128 %======================================================================== 0129 % Projection angles: 0130 % Difference of the k space vector angle and the angles of the ellipses 0131 %======================================================================== 0132 sind = sin(kphi - alpha); 0133 cosd = cos(kphi - alpha); 0134 0135 akphi = sqrt(a.^2 .* cosd.^2 + b.^2 .* sind.^2); 0136 %======================================================================== 0137 % Calculate the signals for each ellips and each k space point kx,ky 0138 % avoid division by zero besselj(1,0) = 0 so set 0/0 = 1 0139 %======================================================================== 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 %======================================================================== 0145 % Sum the signals for all ellipses and reshape the data back to 0146 % the matrix size of the kx and ky input data 0147 %======================================================================== 0148 Sk = sum(Sk,2); 0149 Sk = reshape(Sk, mk,nk); 0150 %======================================================================== 0151 % Show the result, including the reconstruction if default Carthesian 0152 % sampling was used 0153 %======================================================================== 0154 if nargin >= 2 0155 fh = figure; 0156 subplot(1,2,1); 0157 if ~isempty(P) 0158 imagesc(abs(P)); 0159 end; 0160 title('Software head phantom') 0161 axis image 0162 0163 subplot(1,2,2); 0164 imagesc(abs(Sk)); 0165 title('Simulated raw data') 0166 return 0167 else 0168 fh = figure; 0169 subplot(1,3,1); 0170 if ~isempty(P) 0171 imagesc(abs(P)); 0172 end; 0173 axis image 0174 0175 subplot(1,3,2); 0176 imagesc(abs(Sk)); 0177 title('Simulated raw data') 0178 axis image 0179 0180 subplot(1,3,3); 0181 I = fftshift(fft2(Sk)); 0182 imagesc(abs(I)); 0183 axis image 0184 title('2DFFT Reconstructed image') 0185 end; 0186 %======================================================================== 0187 % Suppress output if not asked for any 0188 %======================================================================== 0189 if nargout < 1 0190 Sk = []; 0191 end; 0192 0193 0194 %======================================================================== 0195 % END 0196 %======================================================================== 0197 0198 %======================================================================== 0199 % local function Copied from function phantom.m 0200 % to produce a result without the image toolbox 0201 %======================================================================== 0202 0203 %======================================================================== 0204 function toft=modified_shepp_logan 0205 % 0206 % This head phantom is the same as the Shepp-Logan except 0207 % the intensities are changed to yield higher contrast in 0208 % the image. Taken from Toft, 199-200. 0209 % 0210 % A a b x0 y0 phi 0211 % --------------------------------- 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 %======================================================================== 0223 % end local function toft=modified_shepp_logan 0224 %========================================================================