ctMotionCorrection/testAlgorithms.m at master · ndwork/ctMotionCorrection · 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
function [recon,costs] = testAlgorithms( nDetectors, ...
  detSize, thetas, translations, nCols, nRows, pixSize, varargin)
  % Optional Variables:
  % method:
	%   'GD' for Gradient Descent
  %   'LADMM' for Linearized ADMM
  %   'PC' for Pock Chambolle

  defaultMethod = 'PC';
  defaultInput = 'rand';
  p = inputParser;
  p.addOptional( 'method', defaultMethod );
  p.addOptional('inputMatrix', defaultInput);
  p.parse( varargin{:} );
  method = p.Results.method;
  inputMatrix = p.Results.inputMatrix;

  maxVerticalShift = max(translations(:,1));
  maxHorizontalShift = max(translations(:,2));
  switch inputMatrix

    case 'rand'
      random = rand(100,nRows);
      xTest = rand(nRows,nCols);
      xPad = padImgForRadon( xTest, maxHorizontalShift, maxVerticalShift, ...
        pixSize );

      applyE = @(u) random*u;
      applyET = @(u) random'*u;

      sinogram = applyE(xPad);

      gamma = 0;

    case 'deblur'
      % kernel = fspecial('gaussian',[15 15],7);
      kernel = fspecial('gaussian',[9 9],3); % use this one with the Barbara image below.
      % kernel = randn(15,15);

      kernelTrans = fliplr(flipud(kernel));

      applyA = @(x) imfilter(x,kernel,'circular');  % K is a blur operator.
      applyAT = @(y) imfilter(y,kernelTrans,'circular');

      I = imread('barbara.png');
      I = double(I(:,:,1));
%       I = imresize(I,.5);

      [nRows,nCols] = size(I);

      mn = min(I(:));
      I = I - mn;
      mx = max(I(:));
      I = I/mx;
      xTest = I;

      b = applyA(I);
      noise = (1e-3)*randn(nRows,nCols);
      b = b + noise;
      figure('Name','b')
      imshow(b,[]) % b is a blurry image, ready to be deblurred.
      sinogram = b;

      gamma = 1 / 100000;

      eigValArr_A = eigValArrForCyclicConvOp(kernel,nRows,nCols);
      eigValArr_ATrans = conj(eigValArr_A);

      applyE = @(x) applyCyclicConv2D(x,eigValArr_A);
      applyET = @(x) applyCyclicConv2D(x,eigValArr_ATrans);

  end


  switch method

    case 'GD'     % Gradient Descent
      costs = [];
      recon = ctCorrectForTranslation_GD_test( sinogram, ...
        applyE, applyET, ...
        nDetectors, detSize, thetas, translations, nCols, nRows, ...
        pixSize );

    case 'LADMM'  % Linearized ADMM
      [recon, costs] = ctCorrectForTranslation_LADMM_test( sinogram, ...
        applyE, applyET, ...
        nDetectors, detSize, thetas, translations, nCols, nRows, ...
        pixSize, gamma );

      error = max(norm(recon(:) - xTest(:),2) / norm(xTest(:),2));
      disp(['Error is: ', num2str(error)]);
      disp(['gamma was set to: ', num2str(gamma)])

    case 'PC'     % Pock-Chambolle

      [recon, costs] = ctCorrectForTranslation_PC_test( sinogram, ...
        applyE, applyET, nDetectors, detSize, thetas, translations,...
        nCols, nRows, pixSize, gamma );

      figure; imshow( xTest, [] );
      title('Original Image')

      error = max( norm(recon(:) - xTest(:),2) / norm(xTest(:),2));
      disp(['Error is: ', num2str(error)]);
      disp(['gamma was set to: ', num2str(gamma)])


  end

end