BSReconFramework/main_reconstructProspectivelySubsampled.m at master · IMTtugraz/BSReconFramework · 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
%% add support functions to path
clear all
close all
clc

addpath([pwd,filesep,'supportFunctions']);
addpath([pwd,filesep,'create_mat_file']);
addpath([pwd,filesep,'nabla']);
addpath([pwd,filesep,'tgv']);
addpath([pwd,filesep,'gpu']);
setenv('LD_LIBRARY_PATH','/usr/lib/x86_64-linux-gnu')

%% define parameters
lambda = 70;                    %Regularizaiton parameters
mu = 5e-4;
implementation = 'GPU';         %'GPU': use GPU implementation
                                %'CPU': use CPU implementation
patternType = 'vdrandom';          %specify pattern type
                                %   'full':     fully sampled
                                %   'block':    block pattern in k-space center
                                %   'eliptic':  eliptical pattern in k-space center
                                %   'vdrandom': variable density pattern
                                %   'gauss':    pattern with Gaussian density function

maxitTGV = 1000;                %maximum numbers of iteration
maxitH1  = 1000;


%% load data
FileName = 'gre_BlochSiegert_acc_12x4.dat';
PathName =  './data/';
inFile = [PathName, FileName];
coilSens = 'loadExternal';      	%Coil Sensitivity Maps:
                                    %   'calcFromFullData': coil sensitivities are calculated ouf
                                    %                       of fully sampled data
                                    %   'loadExternal':     coil sensitivities are provided externally
CoilSensPath = './data/smaps_walsh3d_slice.mat';

% format of data3D: [NCol, NLine, NSlice, NPh, NCh]
[raw_data, hdr, PhaseEncDir] = readRAW_BS_3D(inFile);
raw_data = double(raw_data);
disp('Data read completed');

numSlice = hdr.Config.NImagePar;

if strcmp(PhaseEncDir,'LIN')
    dataSize = [hdr.Config.NLinMeas, hdr.Config.NImageCols, hdr.Config.NParMeas];
    imgRes = [hdr.Config.PeFOV, hdr.Config.ReadFoV, hdr.MeasYaps.sSliceArray.asSlice{1}.dThickness]./[dataSize(1:2), numSlice];
elseif strcmp(PhaseEncDir, 'COL')
    dataSize = [hdr.Config.NImageCols, hdr.Config.NLinMeas, hdr.Config.NParMeas];
    imgRes = [hdr.Config.ReadFoV, hdr.Config.PeFOV, hdr.MeasYaps.sSliceArray.asSlice{1}.dThickness]./[dataSize(1:2), numSlice];
end

if max(dataSize ~= [size(raw_data,1), size(raw_data,2), size(raw_data,3)])
    error ('Dimension missmatch')
end


%% Load Parameters
par = loadParameters();
par.impl = implementation;

%% generate pattern

par.pattern = raw_data(:,:,:,:,1) ~= 0;


%% struct par init
[dimY, dimX, dimSlice, Nph, NC] = size(raw_data);
par.NC                = NC;                       % number of coils
par.dimY              = dimY;                     % data dimensions in y
par.dimX              = dimX;                     % data dimensions in x
par.dimSlice          = dimSlice;                 % data dimensions in z

par.y                  = raw_data;
par.PhaseEncDir        = PhaseEncDir;
scaleDim               = max(imgRes(1:2));
par.dx                 = imgRes(1)/scaleDim;
par.dy                 = imgRes(2)/scaleDim;
par.dz                 = imgRes(3)/scaleDim;

par.maxitTGV = maxitTGV;
par.maxitH1  = maxitH1;

par.traj='cart';            %only cartesian implemented
par.coilSens = coilSens;
par.CoilSensPath = CoilSensPath;
if strcmp(par.coilSens, 'calcFromFullData')
    par.CoilSensData = squeeze(raw_data(:,:,:,1,:));
end

clear pattern
disp('Initialization completed');toc


%% starte reco
par.lambda = lambda;
par.mu = mu;

[B1Map, flipAngleMap]= BlochSiegReco(par);

%% remove slice oversampling
sliceDiff = (hdr.Config.NParMeas - numSlice)/2;
flipAngleMap = flipAngleMap(:,:,sliceDiff+1:end-sliceDiff);
B1Map = B1Map(:,:,sliceDiff+1:end-sliceDiff);
pattern = par.pattern(:,:,sliceDiff+1:end-sliceDiff,:,1);

%% write results into data folder
% output: 		flipAngleMap: 	relative flip angle map in % of the nominal flip angle
%               B1Map:          B1+ map in µT
%               pattern:        used subsampling pattern

save([PathName, '/B1Map_prospectivelySubsampled.mat'], 'flipAngleMap','B1Map', 'pattern')