This file is indexed.

/usr/share/octave/packages/image-2.4.1/phantom.m is in octave-image 2.4.1-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
## Copyright (C) 2010 Alex Opie <lx_op@orcon.net.nz>
## Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{P}} = phantom ()
## @deftypefnx {Function File} {@var{P}} = phantom (@var{model})
## @deftypefnx {Function File} {@var{P}} = phantom (@var{E})
## @deftypefnx {Function File} {@var{P}} = phantom (@dots{}, @var{n})
## @deftypefnx {Function File} {[@var{P}, @var{E}]} = phantom (@dots{})
## Create computational phantom head.
##
## A phantom is a known object (either real or purely mathematical) that is
## used for testing image reconstruction algorithms.  The Shepp-Logan phantom
## is a popular mathematical model of a cranial slice, made up of a set of
## overlaying ellipses.  This allows rigorous testing of computed tomography
## (CT) algorithms as it can be analytically transformed with the radon
## transform (see the functions @code{radon} and @code{iradon}).
##
## The phantom @var{P}, is created by overlaying ellipses as defined by the
## matrix @var{E} or one of the standard @var{model}s, in a square of size
## @var{n} by @var{n} (defaults to 256).
##
## The available standard @var{model}s (use the output argument @var{E} to
## inspect the details of the different ellipses) are:
##
## @table @asis
## @item @qcode{"Sheep-Logan"}
## This is the original Sheep-Logan model with 10 ellipses as described in
## Table 1 of @cite{Shepp, Lawrence A., and Benjamin F. Logan. "The Fourier
## reconstruction of a head section." Nuclear Science, IEEE Transactions on
## 21, no. 3 (1974): 21-43.}
##
## @item @qcode{"Modified Shepp-Logan"} (default)
## A modification of the original Shepp-Logan model to give a better contrast,
## as described in Table B.3 of @cite{Toft, Peter Aundal. "The radon
## transform-theory and implementation." PhD diss., Department of Mathematical
## Modelling, Technical University of Denmark, 1996.}
##
## @end table
##
## A 6 column matrix @var{E} can be used to generate a custom image by
## superimposing arbitrary ellipses.  Each row defines a single ellipse, with
## each column for the values of @{I, a, b, x0, y0, phi@}:
##
## @table @abbr
## @item I
##   is the additive intensity of the ellipse
##
## @item a
##    is the length of the major axis
##
## @item b
##    is the length of the minor axis
##
## @item x0
##    is the horizontal offset of the centre of the ellipse
##
## @item y0
##    is the vertical offset of the centre of the ellipse
##
## @item phi
##    is the counterclockwise rotation of the ellipse in degrees,
##    measured as the angle between the x axis and the ellipse major axis.
##
## @end table
##
## The image bounding box in the algorithm is @{[-1, -1], [1, 1]@}, so the
## values of a, b, x0, y0 should all be specified with this in mind.
##
## Example:
##
## @example
## @group
## P = phantom (512);
## imshow (P);
## @end group
## @end example
##
## @seealso{iradon, radon}
## @end deftypefn

function [head, ellipses] = phantom (varargin)

  if (nargin > 2)
    print_usage ()
  endif

  ## Would be really cool if we implemented a 3D phantom as already described
  ## in Cheng Guan Koay, Joelle E. Sarlls, and Evren Ozarslan (2007).
  ## "Three-Dimensional Analytical Magnetic Resonance Imaging Phantom in the
  ## Fourier Domain". Magnetic Resonance in Medicine 58:430 - 436.
  ## The Table 1 on their paper to generate the 3D model, would take 8 columns,
  ## an extra value for z axis coordinates, and extra axis length.
  ## They mention other phantom heads as more canonical 3D head phantoms (read
  ## the introduction)

  ## Defaults
  ellipses  = mod_shepp_logan ();
  n         = 256;

  if (nargin)
    ## Check validity of N
    chk_n = @(x) isnumeric (x) && isscalar (x) && ceil (x) == x;

    in = varargin{1};
    if (ischar (in))
      switch (tolower (in))
        case "shepp-logan",           ellipses = shepp_logan ();
        case "modified shepp-logan",  ellipses = mod_shepp_logan ();
        otherwise
          error ("phantom: unknown MODEL `%s'", in);
      endswitch
    elseif (isnumeric (in) && ndims (in) == 2 && columns (in) == 6)
      ellipses = in;
    elseif (chk_n (in))
      n = in;
      ## If N is the first argument, we can't have more
      if (nargin > 1)
        print_usage ();
      endif
    else
      error ("phantom: first argument must either be MODEL, E, or N");
    endif

    ## If there is a second input argument, must be N
    if (nargin > 1)
      if (chk_n (varargin{2}))
        n = varargin{2};
      else
        error ("phantom: N must be numeric scalar");
      endif
    endif
  endif

  ## Initialize blank image
  head = zeros (n);

  # Create the pixel grid
  xvals = (-1 : 2 / (n - 1) : 1);
  xgrid = repmat (xvals, n, 1);

  for i = 1:rows (ellipses)
    I   = ellipses (i, 1);
    a2  = ellipses (i, 2)^2;
    b2  = ellipses (i, 3)^2;
    x0  = ellipses (i, 4);
    y0  = ellipses (i, 5);
    phi = ellipses (i, 6) * pi / 180;  # Rotation angle in radians

    ## Create the offset x and y values for the grid
    x = xgrid - x0;
    y = rot90 (xgrid) - y0;

    cos_p = cos (phi);
    sin_p = sin (phi);

    ## Find the pixels within the ellipse
    locs = find (((x .* cos_p + y .* sin_p).^2) ./ a2 ...
      + ((y .* cos_p - x .* sin_p).^2) ./ b2 <= 1);

    ## Add the ellipse intensity to those pixels
    head(locs) += I;
  endfor
endfunction

function ellipses = shepp_logan ()
  ## Standard head phantom, taken from Shepp & Logan
  ##
  ## Note that the first element of this matrix, the gray value for the first
  ## ellipse (human skull), has a value of 1.0 even though the paper gives it a
  ## a value of 2.0 (see Table 1 on page 32 and Figure 1 on page 34). This
  ## change is so that the **head** intensity values appear in the range [0 1]
  ## rather than the range [1 2].
  ##
  ##    **The problem with this**
  ##
  ## The background still need an intensity value which is going to be 0. This
  ## means that we can't distinguish between the background and the ventricles
  ## (ellipse "c" and "d" whose intensities are a + b + c and a + b + d, see
  ## Figure 1) since they will have an intensity value of 0 (actually, because
  ## of machine precision the ventricules will be almost 0). But if we didn't
  ## made this change, the ** image** range would be [0 2] with all of the head
  ## details compressed in half of the display range. Also, Matlab seems to be
  ## doing the same.
  persistent ellipses = [ 1     0.69    0.92     0      0         0
                         -0.98  0.6624  0.874    0     -0.0184    0
                         -0.02  0.11    0.31     0.22   0       -18
                         -0.02  0.16    0.41    -0.22   0        18
                          0.01  0.21    0.25     0      0.35      0
                          0.01  0.046   0.046    0      0.1       0
                          0.01  0.046   0.046    0     -0.1       0
                          0.01  0.046   0.023   -0.08  -0.605     0
                          0.01  0.023   0.023    0     -0.606     0
                          0.01  0.023   0.046    0.06  -0.605     0];
endfunction

function ellipses = mod_shepp_logan ()
  ## Modified version of Shepp & Logan's head phantom, adjusted to improve
  ## contrast.  Taken from Peter Toft PhD thesis, Table B.3
  persistent ellipses = [ 1.0   0.69    0.92     0.0    0.0       0
                         -0.8   0.6624  0.874    0.0   -0.0184    0
                         -0.2   0.11    0.31     0.22   0.0     -18
                         -0.2   0.16    0.41    -0.22   0.0      18
                          0.1   0.21    0.25     0.0    0.35      0
                          0.1   0.046   0.046    0.0    0.1       0
                          0.1   0.046   0.046    0.0   -0.1       0
                          0.1   0.046   0.023   -0.08  -0.605     0
                          0.1   0.023   0.023    0.0   -0.606     0
                          0.1   0.023   0.046    0.06  -0.605     0];
endfunction

%!demo
%! P = phantom (512);
%! imshow (P);