function SI2con_SRFLASH_roi(data1, info, bw, skip, PhaseEncSteps, TR_real)
% test conversion to a value linear w.r.t. concentration on a ROI
% input args:
%   data1, data2, info - input data
%   bw - mask of the testing roi
%   skip - how many frames to skip
%   PhaseEncSteps - number of phase encoding steps before the acquisition of the middle
%     k-space line
%   TR_real - TR of the acquisition sequence in [ms]




%% precontrast sequences - extract imagesc for estimation of M0 and E10

seqs = info.acq.kM0Range(2)-info.acq.kM0Range(1)+1;   % number of precontrast seqs.
slices = info.sections;

tis = 1:3000; % TI axis

%% get image size, and number of frames
frames = 0;
for s=1:seqs
    f = length(data1{info.acq.kM0Range(1)+s-1});
    if f>frames
        frames=f;
    end
end
X = size(data1{1}{1},2);
Y = size(data1{1}{1},1);

%% ROI
numpixels = sum(sum(bw));


si = zeros(1,50);
figure(100);
hold off;
plot(tis(1),si(1));
hold on;

%img = cell(seqs,frames);   % sequence , frame
% img = zeros(Y,X,seqs,frames);   % sequence , frame
% SI =  zeros(seqs,frames);


%%
numframes = 0; % number of all precontrast images
xdata = zeros(0);
for seq = info.acq.kM0Range(1):info.acq.kM0Range(2) %1:seqs
    frn = length(data1{seq}) - skip; % number of frames of the given sequence - skipframes
    xdata = [xdata double((info.acq.TI(seq) - PhaseEncSteps*TR_real)*ones(1,frn)) ];
    numframes = numframes + frn;
end

img = zeros(Y,X,numframes);   % sequence , frame
% SI = zeros(numframes);
ind = 1;
FA = info.acq.FA(info.acq.kM0Range(1));
for seq = info.acq.kM0Range(1):info.acq.kM0Range(2) %1:seqs
    frn = length(data1{seq}) - skip; % number of frames of the given sequence - skipframes
    for frame=1:frn
        img(:,:,ind) = double( data1{seq}{frame+skip} );
        
        ti = xdata(ind);
        si = sum(sum(double(bw).* img(:,:,ind))) / numpixels;
        ydata(ind)=si;
        %figure(2);imagesc(im + im.*double(bw)*1); colormap(gray);
        %title(sprintf([precontrast_seq_mat{seq} '   frame: %d'],slice,frame));
        figure(100);
        switch frame
            case 1
                plot(ti,si,'.');
            case 2
                plot(ti,si,'*');
            case 3
                plot(ti,si,'x');
            otherwise
                plot(ti,si,'o');
        end
        hold on
        ind = ind + 1;
    end
    if info.acq.FA(seq) ~= FA
        disp('Error: different FA values in precontrast acquisition! ')
        return
    end
end


%% curve fitting to estimate E10 and M0 for ROI:

% range = info.acq.kM0Range(1):info.acq.kM0Range(2);
%xdata = double(repmat(info.acq.TI(range)', [1,frames-info.acq.skip_frames]));
% ydata = SI(:,info.acq.skip_frames+1:end);
nm1 = PhaseEncSteps;    % n-1

% get initial value for M0:
M0_odhad = max(ydata);

% curve fitting:
% [T10 M0] = SI2con_fit_SRFLASH(xdata, ydata, TR_real, FA, nm1, M0_odhad);
[T10 M0] = SI2con_fit_SRFLASH(xdata, ydata, TR_real, FA, nm1);
%T10 = -info.acq.TR/log(E10);

%% construct the fitted curve:
E10 = exp(-tis/T10);
%y = M0 * (1-E10);

TR = TR_real; %info.acq.TR_real;
% FA = info.acq.FA(range(1));

E = exp(-TR/T10);
a = cos(FA/180*pi) .* E;

y_auto = M0 * ( (1-E10).*(a.^nm1) + (1-E).*(1-a.^nm1)./(1-a));

%figure(1);plot(xdata, ydata, 'o', tis,y,'-');
%figure(1);
%plot(tis,y,'-');
% title(sprintf('fitting results: M0: %f,   T10: %f ms',M0, T10))

M0a=M0;
T10a=T10;

%%  Complete search optimization:

disp('manual optimisation res:');

T1_weight = 1000;

M0_a = M0_odhad*0.1:M0_odhad/10:M0_odhad*50;  %1000:100:6000; %M0-3000:100:M0+5000; %M0*(0.1:0.1:10);
%E10_a = (E10-0.01:0.001:1)*E1_weight; %E10*(0.1:0.1:10)*E1_weight;
T10_a = (50:10:2000) * T1_weight;
Mc = length(M0_a);
%Ec = length(E10_a);
Ec = length(T10_a);
SI_mat = zeros(Mc,Ec);
for m=1:Mc
    for e=1:Ec
        %y=SI2con_precontrast([M0_a(m) E10_a(e)],xdata);
        %y=SI2con_precontrast_T1_v2([M0_a(m) T10_a(e)],xdata,info.acq.TR, E1_weight);
        y= SI2con_precontrast_SRFLASH_eval([M0_a(m) T10_a(e)],xdata,T1_weight,TR,FA,nm1);
        SI_mat(m,e) = sum(sum((y-ydata).^2));
    end
end
% figure(33);
% surf(T10_a/T1_weight, M0_a, SI_mat);
% ylabel('M0');
% xlabel('T10');

[a, inds] =min(SI_mat); %minima along M0 (columns)
[a, ind] =min(a); %minima along T10
M0 = M0_a(inds(ind))
T10 = T10_a(ind)/T1_weight

disp(sprintf('complete search results: M0: %f,   T10: %f ms',M0, T10))

E10 = exp(-tis/T10);
%y = M0 * (1-E10);

TR = TR_real;
% FA = info.acq.FA(range(1));

E = exp(-TR/T10);
a = cos(FA/180*pi) .* E;

y_man = M0 * ( (1-E10).*(a.^nm1) + (1-E).*(1-a.^nm1)./(1-a));

figure(100);plot(tis,y_auto,'b-',   tis,y_man,'g:');
% figure(100);plot(tis,y_auto,'-');
legend('Newton optimization', '"Complete search" optimization');
title(sprintf('automatic fitting results: M0: %f,   T10: %f ms',M0a, T10a))

hold off

% figure(6); hold off


end




