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
|
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 ASTRA Toolbox
%
% Copyright: 2010-2021, imec Vision Lab, University of Antwerp
% 2014-2021, CWI, Amsterdam
% License: Open Source under GPLv3
% Contact: astra@astra-toolbox.com
% Website: http://www.astra-toolbox.com/
%------------------------------------------------------------------------
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);
|
