summaryrefslogtreecommitdiffstats
path: root/matlab/tools/rebin_fan2par.m
blob: da20932e50801522e6f403a56d1eb292644a42e9 (plain)
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
function F = rebin_fan2par(RadonData, BetaDeg, D, thetaDeg)

%------------------------------------------------------------------------
% F = rebin_fan2par(RadonData, BetaDeg, D, thetaDeg)
%
% Deze functie zet fan beam data om naar parallelle data, door interpolatie
% (fast resorting algorithm, zie Kak en Slaney)
% Radondata zoals altijd: eerste coord gamma , de rijen
%                          tweede coord beta, de kolommen, beide hoeken in
%                          radialen
% PixPProj: aantal pixels per projectie (voor skyscan data typisch 1000)
% BetaDeg: vector met alle draaihoeken in graden
% D: afstand bron - rotatiecentrum in pixels, dus afstand
% bron-rotatiecentrum(um) gedeeld door image pixel size (um).
% thetaDeg: vector met gewenste sinogramwaarden voor theta in graden
%       de range van thetaDeg moet steeds kleiner zijn dan die van betadeg
% D,gamma,beta, theta zoals gebruikt in Kak & Slaney
%------------------------------------------------------------------------
%------------------------------------------------------------------------
% This file is part of the
% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
%
% Copyright: iMinds-Vision Lab, University of Antwerp
% License: Open Source under GPLv3
% Contact: mailto:astra@ua.ac.be
% Website: http://astra.ua.ac.be
%------------------------------------------------------------------------
% $Id$

NpixPProj = size(RadonData,1);  % aantal pixels per projectie
%if mod(size(Radondata,1),2)==0
%    NpixPProjNew=NpixPProj+1;
%else
    NpixPProjNew = NpixPProj;
%end

%% FAN-BEAM RAYS

% flip sinogram, why?
RadonData = flipdim(RadonData,2);  %  matlab gebruikt tegengestelde draairichting (denkik) als skyscan, of er is een of andere flipdim geweest die gecorrigeerd moet worden))

% DetPixPos: distance of each detector to the ray through the origin (theta)
DetPixPos = -(NpixPProj-1)/2:(NpixPProj-1)/2;  % posities detectorpixels

% GammaStralen: alpha's? (result in radians!!)
GammaStralen = atan(DetPixPos/D); % alle met de detectorpixelposities overeenkomstige gammahoeken

% put beta (theta) and gamma (alpha) for each ray in 2D matrices
[beta gamma] = meshgrid(BetaDeg,GammaStralen);

% t: minimal distance between each ray and the ray through the origin
t = D*sin(gamma); % t-waarden overeenkomstig met de verschillende gamma's

theta = gamma*180/pi + beta;  % theta-waarden in graden overeenkomstig met verschillende gamma en beta waarden

%% PARALLEL BEAM RAYS

% DetPixPos: distance of each detector to the ray through the origin (theta)
DetPixPos = -(NpixPProjNew-1)/2:(NpixPProjNew-1)/2;  % posities detectorpixels

% GammaStralen: alpha's? (result in radians!!)
GammaStralenNew = atan(DetPixPos/D); % alle met de detectorpixelposities overeenkomstige gammahoeken

% put beta (theta) and gamma (alpha) for each ray in 2D matrices
[~, gamma] = meshgrid(BetaDeg,GammaStralenNew);

% t: minimal distance between each ray and the ray through the origin
tnew = D * sin(gamma); % t-waarden overeenkomstig met de verschillende gamma's

% calculate new t
step = (max(tnew)-min(tnew)) / (NpixPProjNew-1);
t_para = min(tnew):step:max(tnew);

[thetaNewCoord tNewCoord] = meshgrid(thetaDeg, t_para);

%% Interpolate
Interpolant = TriScatteredInterp(theta(:), t(:), RadonData(:),'nearest');
F = Interpolant(thetaNewCoord,tNewCoord);