phantom3d.m 8.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234
  1. % PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom
  2. % P = PHANTOM3D(DEF,N) generates a 3D head phantom that can
  3. % be used to test 3-D reconstruction algorithms.
  4. %
  5. % DEF is a string that specifies the type of head phantom to generate.
  6. % Valid values are:
  7. %
  8. % 'Shepp-Logan' A test image used widely by researchers in
  9. % tomography
  10. % 'Modified Shepp-Logan' (default) A variant of the Shepp-Logan phantom
  11. % in which the contrast is improved for better
  12. % visual perception.
  13. %
  14. % N is a scalar that specifies the grid size of P.
  15. % If you omit the argument, N defaults to 64.
  16. %
  17. % P = PHANTOM3D(E,N) generates a user-defined phantom, where each row
  18. % of the matrix E specifies an ellipsoid in the image. E has ten columns,
  19. % with each column containing a different parameter for the ellipsoids:
  20. %
  21. % Column 1: A the additive intensity value of the ellipsoid
  22. % Column 2: a the length of the x semi-axis of the ellipsoid
  23. % Column 3: b the length of the y semi-axis of the ellipsoid
  24. % Column 4: c the length of the z semi-axis of the ellipsoid
  25. % Column 5: x0 the x-coordinate of the center of the ellipsoid
  26. % Column 6: y0 the y-coordinate of the center of the ellipsoid
  27. % Column 7: z0 the z-coordinate of the center of the ellipsoid
  28. % Column 8: phi phi Euler angle (in degrees) (rotation about z-axis)
  29. % Column 9: theta theta Euler angle (in degrees) (rotation about x-axis)
  30. % Column 10: psi psi Euler angle (in degrees) (rotation about z-axis)
  31. %
  32. % For purposes of generating the phantom, the domains for the x-, y-, and
  33. % z-axes span [-1,1]. Columns 2 through 7 must be specified in terms
  34. % of this range.
  35. %
  36. % [P,E] = PHANTOM3D(...) returns the matrix E used to generate the phantom.
  37. %
  38. % Class Support
  39. % -------------
  40. % All inputs must be of class double. All outputs are of class double.
  41. %
  42. % Remarks
  43. % -------
  44. % For any given voxel in the output image, the voxel's value is equal to the
  45. % sum of the additive intensity values of all ellipsoids that the voxel is a
  46. % part of. If a voxel is not part of any ellipsoid, its value is 0.
  47. %
  48. % The additive intensity value A for an ellipsoid can be positive or negative;
  49. % if it is negative, the ellipsoid will be darker than the surrounding pixels.
  50. % Note that, depending on the values of A, some voxels may have values outside
  51. % the range [0,1].
  52. %
  53. % Example
  54. % -------
  55. % ph = phantom3d(128);
  56. % figure, imshow(squeeze(ph(64,:,:)))
  57. %
  58. % Copyright 2005 Matthias Christian Schabel (matthias @ stanfordalumni . org)
  59. % University of Utah Department of Radiology
  60. % Utah Center for Advanced Imaging Research
  61. % 729 Arapeen Drive
  62. % Salt Lake City, UT 84108-1218
  63. %
  64. % This code is released under the Gnu Public License (GPL). For more information,
  65. % see : http://www.gnu.org/copyleft/gpl.html
  66. %
  67. % Portions of this code are based on phantom.m, copyrighted by the Mathworks
  68. %
  69. % =========================================================================
  70. %% Phantom Code
  71. function [p,ellipse]=phantom3d(varargin)
  72. [ellipse,n] = parse_inputs(varargin{:});
  73. p = ones([n n n]) .* 0.0357776; % Attenuation coefficient of PMMA (1/mm)
  74. %p = zeros([n n n]);
  75. rng = ((0:n-1)-(n-1)/2) / ((n-1)/2);
  76. [x,y,z] = meshgrid(rng,rng,rng); % Coordinate system of a phantom; starting point at the phantom's centre
  77. coord = [flatten(x); flatten(y); flatten(z)];
  78. p = flatten(p);
  79. for k = 1:size(ellipse,1)
  80. A = ellipse(k,1); % Amplitude change for this ellipsoid
  81. asq = ellipse(k,2)^2; % a^2
  82. bsq = ellipse(k,3)^2; % b^2
  83. csq = ellipse(k,4)^2; % c^2
  84. x0 = ellipse(k,5); % x offset
  85. y0 = ellipse(k,6); % y offset
  86. z0 = ellipse(k,7); % z offset
  87. phi = ellipse(k,8)*pi/180; % first Euler angle in radians
  88. theta = ellipse(k,9)*pi/180; % second Euler angle in radians
  89. psi = ellipse(k,10)*pi/180; % third Euler angle in radians
  90. cphi = cos(phi);
  91. sphi = sin(phi);
  92. ctheta = cos(theta);
  93. stheta = sin(theta);
  94. cpsi = cos(psi);
  95. spsi = sin(psi);
  96. % Euler rotation matrix
  97. alpha = [cpsi*cphi-ctheta*sphi*spsi cpsi*sphi+ctheta*cphi*spsi spsi*stheta;
  98. -spsi*cphi-ctheta*sphi*cpsi -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta;
  99. stheta*sphi -stheta*cphi ctheta];
  100. % Rotated ellipsoid coordinates
  101. coordp = alpha*coord;
  102. idx = find((coordp(1,:)-x0).^2./asq + (coordp(2,:)-y0).^2./bsq + (coordp(3,:)-z0).^2./csq <= 1);
  103. p(idx) = p(idx) + A;
  104. end
  105. p = reshape(p,[n n n]);
  106. return;
  107. function out = flatten(in)
  108. out = reshape(in,[1 prod(size(in))]);
  109. return;
  110. function [e,n] = parse_inputs(varargin)
  111. % e is the m-by-10 array which defines ellipsoids
  112. % n is the size of the phantom brain image
  113. n = 128; % The default size
  114. e = [];
  115. defaults = {'shepp-logan', 'modified shepp-logan', 'yu-ye-wang', 'aleks'};
  116. for i=1:nargin
  117. if ischar(varargin{i}) % Look for a default phantom
  118. def = lower(varargin{i});
  119. idx = strmatch(def, defaults);
  120. if isempty(idx)
  121. eid = sprintf('Images:%s:unknownPhantom',mfilename);
  122. msg = 'Unknown default phantom selected.';
  123. error(eid,'%s',msg);
  124. end
  125. switch defaults{idx}
  126. case 'shepp-logan'
  127. e = shepp_logan;
  128. case 'modified shepp-logan'
  129. e = modified_shepp_logan;
  130. case 'yu-ye-wang'
  131. e = yu_ye_wang;
  132. case 'aleks'
  133. e = aleks;
  134. end
  135. elseif numel(varargin{i})==1
  136. n = varargin{i}; % a scalar is the image size
  137. elseif ndims(varargin{i})==2 && size(varargin{i},2)==10
  138. e = varargin{i}; % user specified phantom
  139. else
  140. eid = sprintf('Images:%s:invalidInputArgs',mfilename);
  141. msg = 'Invalid input arguments.';
  142. error(eid,'%s',msg);
  143. end
  144. end
  145. % ellipse is not yet defined
  146. if isempty(e)
  147. e = modified_shepp_logan;
  148. end
  149. return;
  150. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  151. % Default head phantoms: %
  152. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  153. function e = shepp_logan
  154. e = modified_shepp_logan;
  155. e(:,1) = [1 -.98 -.02 -.02 .01 .01 .01 .01 .01 .01];
  156. return;
  157. function e = aleks
  158. e = readmatrix("Mikrokalcinacije.txt")
  159. return;
  160. function e = modified_shepp_logan
  161. %
  162. % This head phantom is the same as the Shepp-Logan except
  163. % the intensities are changed to yield higher contrast in
  164. % the image. Taken from Toft, 199-200.
  165. %
  166. % A a b c x0 y0 z0 phi theta psi
  167. % -----------------------------------------------------------------
  168. e = [ 1 .6900 .920 .810 0 0 0 0 0 0
  169. -.8 .6624 .874 .780 0 -.0184 0 0 0 0
  170. -.2 .1100 .310 .220 .22 0 0 -18 0 10
  171. -.2 .1600 .410 .280 -.22 0 0 18 0 10
  172. .1 .2100 .250 .410 0 .35 -.15 0 0 0
  173. .1 .0460 .046 .050 0 .1 .25 0 0 0
  174. .1 .0460 .046 .050 0 -.1 .25 0 0 0
  175. .1 .0460 .023 .050 -.08 -.605 0 0 0 0
  176. .1 .0230 .023 .020 0 -.606 0 0 0 0
  177. .1 .0230 .046 .020 .06 -.605 0 0 0 0 ];
  178. return;
  179. function e = yu_ye_wang
  180. %
  181. % Yu H, Ye Y, Wang G, Katsevich-Type Algorithms for Variable Radius Spiral Cone-Beam CT
  182. %
  183. % A a b c x0 y0 z0 phi theta psi
  184. % -----------------------------------------------------------------
  185. e = [ 1 .6900 .920 .900 0 0 0 0 0 0
  186. -.8 .6624 .874 .880 0 0 0 0 0 0
  187. -.2 .4100 .160 .210 -.22 0 -.25 108 0 0
  188. -.2 .3100 .110 .220 .22 0 -.25 72 0 0
  189. .2 .2100 .250 .500 0 .35 -.25 0 0 0
  190. .2 .0460 .046 .046 0 .1 -.25 0 0 0
  191. .1 .0460 .023 .020 -.08 -.65 -.25 0 0 0
  192. .1 .0460 .023 .020 .06 -.65 -.25 90 0 0
  193. .2 .0560 .040 .100 .06 -.105 .625 90 0 0
  194. -.2 .0560 .056 .100 0 .100 .625 0 0 0 ];
  195. return;