function [R2map, M0map, ErrorMap] = mappT2(CPMGvolume, CPMGinfo)


R2map = zeros(size(CPMGvolume(:,:,1,1)),'single');
M0map = zeros(size(CPMGvolume(:,:,1,1)),'single');
ErrorMap = zeros(size(CPMGvolume(:,:,1,1)),'single');

R2vec = 1:0.1:40;

%% Inital fit
R2slice = zeros([size(CPMGvolume(:,:,1,1)) length(R2vec)],'single');

OneVolume = ones(size(CPMGvolume(:,:,1,1)),'single');

fitIDX = 0;

for R2 = R2vec
    fitIDX = fitIDX + 1;
    R2slice(:,:,fitIDX) = R2*OneVolume;
end

estCPMG = modelCPMG(R2slice,CPMGinfo);
c = mean(CPMGvolume./estCPMG,4);
c(c<200) = 200;
c(c>5000) = 5000;

residual = sum(sqrt((c.*estCPMG - CPMGvolume).^2),4);


[~, minIDX] = min(residual,[],3,'linear');

R2map(:,:,1) = R2slice(minIDX);
ErrorMap(:,:,1) = residual(minIDX);
M0map(:,:,1) = c(minIDX);

%% fine tune


fitIDX = 0;

R2vecAdd = -0.2:0.01:0.2;

for R2add = R2vecAdd
    fitIDX = fitIDX + 1;
    R2slice(:,:,fitIDX) = R2map+R2add;
end

estCPMG = modelCPMG(R2slice,CPMGinfo);
c = mean(CPMGvolume./estCPMG,4);
c(c<200) = 200;
c(c>5000) = 5000;

residual = sum(sqrt((c.*estCPMG - CPMGvolume).^2),4);


[~, minIDX] = min(residual,[],3,'linear');


R2map(:,:,1) = R2slice(minIDX);
ErrorMap(:,:,1) = residual(minIDX);
M0map(:,:,1) = c(minIDX);


R2map(:,:,1) = R2slice(minIDX);
R2map(M0map==200) = 0;

ErrorMap(:,:,1) = residual(minIDX);
M0map(:,:,1) = c(minIDX);

end

function estCPMG = modelCPMG(R2slice,info)
% Sub function for modelling CPMG signal


TEvec = info.TE;

estCPMG = zeros([size(R2slice) length(TEvec)]);

for TEidx = 1:length(TEvec)

    estCPMG(:,:,:,TEidx) = exp(-R2slice.*TEvec(TEidx));

end
end