ctMotionCorrection/run_ctCorrectForRotAndTrans.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

function run_ctCorrectForRotAndTrans
  clear; close all;
  addpath(genpath('.'));

  % Reconstruction parameters
  method = 'PC';    % Options: GD, PC, LADMM
  cy = 0;   nRows=32;
  cx = 0;   nCols=32;
  pixSize = 0.001; % meters / pixel

  datacase = 1;
  switch datacase
    case 1
      im = phantom();

    case 2
      im = double( imread( 'lena.png' ) );
  end
  im = imresize( im, [nCols nRows], 'bilinear' );

  figure; imshow( imresize(im,10,'nearest'), [] );
  title('original');  drawnow;

  detSize = 0.001;
  dTheta = 1 * pi/180;
  thetas = 0:dTheta:pi-dTheta;
  nThetas = numel(thetas);

  nDetectors = nCols*2;

  nonzeroTranslations = 0;
  if nonzeroTranslations
    maxVerticalShift = 0.01; % in meters
    maxHorizontalShift = 0.02; % in meters
    translations = zeros( nThetas, 2 );
    translations(:,1) = linspace(0,maxVerticalShift,nThetas);
    translations(:,2) = linspace(0,maxHorizontalShift,nThetas);
  else
    maxVerticalShift = 0; % in meters
    maxHorizontalShift = 0; % in meters
    translations = zeros( nThetas, 2 );
  end

  nonzeroRotations = 1;
  if nonzeroRotations
    minRotation = -10 * pi/180;
    maxRotation = 10 * pi/180;
    rotations = linspace(minRotation,maxRotation,nThetas);
  end

  im = padImgForRadon( im, maxHorizontalShift, maxVerticalShift, ...
      pixSize );
  [nRows,nCols] = size(im);

  sinogram = radonWithRotAndTrans( im, pixSize, nDetectors, detSize, ...
   thetas, rotations, translations );

  profile on
  tic;
  [recon,costs] = ctCorrectForRotAndTrans( sinogram, nDetectors, ...
    detSize, thetas, rotations, translations, nCols, nRows, pixSize, ...
    'method', method );
  timeTaken = toc;
  profile off
  profile viewer

  disp(['Time taken: ', num2str(timeTaken)]);
  figure; imshow( imresize(recon,10,'nearest'), [] );
  title('Reconstructed image');

  figure; plot( costs, 'LineWidth', 2 );
  xlabel('Iteration'); ylabel('Cost Function');
end