ctCodes/ctRadon.m at master · ndwork/ctCodes · GitHub

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

function sinogram = ctRadon( img, delta, nDetectors, dSize, thetas, ...
  varargin )
  % img:  2D array - will take the Radon transform of this image
  % delta: horizontal and vertical size of pixel (assumed square)
  % nDetectors: the number of detectors
  % thetas: a 1D array, each element is the angle that corresponds to row
  %    radon domain
  % type (optional): 'iso' or 'fast'
  %    'iso' (default) uses a rotation that's an isometry
  %    'fast' faster implementation

  defaultType = 'fast';
  expectedTypes = { 'iso', 'fast' };

  p = inputParser;
  p.addRequired('img', @(x) ndims(x)==2);
  p.addRequired('delta',@isnumeric);
  p.addRequired('nDetectors',@isnumeric);
  p.addRequired('dSize',@isnumeric);
  p.addRequired('thetas');
  p.addOptional('type',defaultType, ...
    @(x) any(validatestring(x,expectedTypes)) );

  p.parse( img, delta, nDetectors, dSize, thetas, ...
    varargin{:} );
  inputs = p.Results;
  type = inputs.type;

  if strcmp( type, 'fast' )
    sinogram = ctRadonFast( img, delta, nDetectors, dSize, thetas );
  elseif strcmp( type, 'iso' )
    sinogram = ctRadonIso( img, delta, nDetectors, dSize, thetas );
  end
end


function sinogram = ctRadonFast( img, delta, nDetectors, dSize, thetas )
  dOffset=0;  % center channel offset

  nTheta = numel(thetas);

  dLocs = ( [0:nDetectors-1] - floor(0.5*nDetectors) ) * dSize - dOffset;

  thetas_deg = thetas * 180/pi;

  Ny = size( img, 1 );  halfY = Ny/2;
  Nx = size( img, 2 );  halfX = Nx/2;
  xs = ones(Ny,1) * (1:Nx);
  ys = (1:Ny)' * ones(1,Nx);
  xs = xs - halfX;
  ys = ys - halfY;
  radiusMask = sqrt( xs.*xs + ys.*ys ) < min(Nx/2,Ny/2);
  radiusImg = img .* radiusMask;

  if mod( Nx, 2 )==0
    locs = ( ([0:Nx-1]) - 0.5*Nx + 0.5 ) * delta;
  else
    locs = ( ([0:Nx-1]) - floor(0.5*Nx) ) * delta;
  end

  sinogram = zeros( nTheta, nDetectors );
  parfor th=1:numel(thetas)
    theta = thetas_deg(th);
    rotImg = imrotate( radiusImg, theta, 'bilinear','crop' );
    sumResult = sum( rotImg, 1 ) * delta;

    interped = interp1( locs, sumResult, dLocs,'linear',0 );

    sinogram(th,:) = interped;
    if mod(th,10)==0 disp(['ctRadon Theta: ', num2str(th), ' of ', ...
        num2str(numel(thetas)) ]); end;
  end

end

function [sinogram B] = ctRadonIso( img, delta, nDetectors, dSize, thetas )
  dOffset=0;  % center channel offset

  nTheta = numel(thetas);

  dLocs = ( [0:nDetectors-1] - floor(0.5*nDetectors) ) * dSize - dOffset;

  Ny = size( img, 1 );  halfY = Ny/2;
  Nx = size( img, 2 );  halfX = Nx/2;
  xs = ones(Ny,1) * (1:Nx);
  ys = (1:Ny)' * ones(1,Nx);
  xs = xs - halfX;
  ys = ys - halfY;
  radiusMask = sqrt( xs.*xs + ys.*ys ) < min(Nx/2,Ny/2);
  radiusImg = img .* radiusMask;

  if mod( Nx, 2 )==0
    locs = ( ([0:Nx-1]) - 0.5*Nx + 0.5 ) * delta;
  else
    locs = ( ([0:Nx-1]) - floor(0.5*Nx) ) * delta;
  end

  % create the interpolation matrix
  B = zeros(nDetectors, Nx);
  for i=1:Nx
    tmp = zeros(1,Nx);
    tmp(i) = 1;
    interped = interp1( locs, tmp, dLocs,'linear',0 );
    B(:,i) = interped;
  end


  sinogram = zeros( nTheta, nDetectors );
  parfor th=1:numel(thetas)
    theta = thetas(th);
    rotImg = isoRot( radiusImg, theta );
    sumResult = sum( rotImg, 1 ) * delta;

    %interped = interp1( locs, sumResult, dLocs,'linear',0 );
    interped = B * sumResult';

    sinogram(th,:) = interped;
    if mod(th,10)==0 disp(['ctRadon Theta: ', num2str(th), ' of ', ...
        num2str(numel(thetas)) ]); end;
  end

end