[Iout,PHt]=perfusion_ph(sigma,Nt,tsc,vp,Kt,Kep,varargin )

function [Iout,PHt]=perfusion_ph(sigma,Nt,tsc,vp,Kt,Kep,varargin )

%PERFUSION_PH
%
% Heart perfussion MRI phantom creation with one scar
%
% OUTPUT:
%
%	Iout: perfusion secuence
%       PHt:   original segmented image 
%
% INPUT:
%
%	sigma:	std of noise (in each coil)
%		if sigma is a matrix is the covariance matrix for multiple coils
%               (correlation betwwen coils assumed)
%	Nt:	number of scans in time (>12)
%	tsc:	time between scans (in seconds)
%	vp:	valvular volume ratio ([0-1]) [default 0]
%	Kt:	transfer constant
%	Kep:	Inverse transfer constant
%
% vp, Kt, Kep can be escalar (for LV tissue) or vector:
%
%	[LV_tissue scar]
%
%
% OPTIONAL INPUT:
%
%	number of coils: 	'singlecoil'
%				'multiplecoil' + [Ncoils] number of coils
%
%	parallel imaging:       (only if multiple coils)
%				'sense'     SENSE simulation 
%                                           (NOT IMPLEMENTED IN THIS VERSION)
%				'grappa'    GRAPPA simulation
%                                           (NOT IMPLEMENTED IN THIS VERSION)
%
%	correlation of noise: 	'corr' (spatially correlated) 
%				'No corr' (not correlated)
%
%	Partial Volume Effect	'PVE' + [sigma] (variance of blur window)
%				'No PVE'
%
%	Sensitivity field	'field'
%
%	
%	Temporal Movement:	'tmove'+[Val] (it will need further registration)
%					Val (float): level of motion: [1-5];
%
%	Maximum Value 		'xmax'+[val] (Default 255)
%	AIF simulation:		'AIF_sim' (use mathematical model)
%				'AIF_file' (use real data from file)
%
% DEFAULT: Rician, correlated, PVE (1.5), AIF_SIM
%
% USAGE:
%
%	[Iout PHt]=perfusion_ph(8,25,1.8,0,[0.5,0.5],[0.3 0.4],'No PVE','tmove',1);
%
% Barcelona, 15/04/2010
% Last modification: 09/02/2011
% Santiago Aja Fernandez
% www.lpi.tel.uva.es/~santi

%PARSE INPUTS------------------------------
%[ncorr,blur,field,sense,tmove, xmax, AIF_SIM]
addpath sources
[Ncoils,ncorr,blur,field,parallel,tmove,xmax,AIF_SIM] = parse_inputs(varargin{:});

%TISSUES----------------------------------------

% AREAS in PHANTOM (indices):
% 	Blood (LV):		255
%       Blood (RV):		247
%	"Around blood": 	128
%	Left Ventricle wall:  	38
%	Scar:		 	26
%	R. ventricle:	 	65
%	Air (lung):	  	0
%	Inner space:	 	13
%	Other blood:		240
%	Other tissues:  	204

PH=imread('ph_heart_1.png');


%IF SENSE length and width must be even
if (parallel==1)||(parallel==2)
   [PHx,PHy]=size(PH);
   if rem(PHx,2)
	PH=PH(1:(end-1),:);
   end
   if rem(PHy,2)
	PH=PH(:,1:(end-1));
   end
end

Blood=(PH==255);
Blood2=(PH==247);
Blood3=(PH==128);
Blood4=(PH==240);
LV=(PH==38);
Scar=(PH==26);
RV=(PH==65);
Lung=(PH==0);
Inter=(PH==13);
Other=(PH==204);

%Original Levels (Before perfusion)

S0=zeros(size(PH));
S0(Blood)=32;
S0(Blood2)=32;
S0(Blood3)=32;
S0(LV)=32;
S0(RV)=32;
S0(Scar)=32;
S0(Inter)=28;
S0(Other)=150;
S0(Blood4)=240;

%M0 Magnetization constant
M0=255;

%if xmax>0
%   S0=S0./255.*xmax;
%end


%============================================================
%IMAGE S0 (Before Perfussion)================================
%============================================================

%NOISE+DISTORTION--------------------------------------------

%Kind of noise:
if Ncoils==1
   Knoise=1;   %Rician
elseif Ncoils>1
   Knoise=parallel+2;
   %Select kind of noise
end

[Mx My]=size(PH);


S0n=distort_sg(S0,sigma,blur,field, Knoise,ncorr,Ncoils);

%============================================================
%PERFUSION IMAGES================================
%============================================================

%Perfusion Params-----------
%if Nt<=10
%	error('Nt>10');
%end


Falldown=0.2; % -->  % of maximum value of decay
delay_LV_RV=0; % -->  Delay of LV and RV in blood filling
delay_LV_RV=round(delay_LV_RV./tsc);
delay_perf=5;
delay_perf=round(delay_perf./tsc)+delay_LV_RV;

t=tsc:tsc:tsc*(Nt);
r1=4.5;
TR=2.7290;
%TR=3;
%FlipAngle= pi/4; %45º
% Simulate c_a(t) --> AIF--------------------
%
%Initial T1 value: 1./T10
%REAL VALUE
T10=-1./TR.*log(1-S0./M0);
%ESTIMATED VALUE
T10e=-1./TR.*log(1-S0n./M0);

%===========================================
% AIF------------------------------
%===========================================

%IF SIMULATED------------------------
if AIF_SIM

F1=10;
F2=1.4;
a1=3;
b1=1;
a2=7;
b2=4;

%cat_RV=F1.*gampdf(t,a1,b1)+F2.*gampdf(t,a2,b2);
%Gamma Model (one gamma)
cat_RV=F1.*gampdf(t,a1,b1);
cat_LV=[zeros([1,delay_LV_RV]) cat_RV];
%LOADED from file
else %IF AIF_FILE
        error('Not implemented yet');
	load s_blood;
	% I_RV: intensity right ventricle
	% I_LV: intensity left ventricle
	[MRV T1RV]=T1_estimate(I_RV,FlipAngle,TR,1);
        [MLV T1LV]=T1_estimate(I_LV,FlipAngle,TR,1);
	cat_RV=((1./T1RV)-(1/T1RV(1)))/r1;
	cat_LV=((1./T1LV)-(1/T1LV(1)))/r1;
end

%Simulate c(t) 

if length(Kep)==1
	Kep=[Kep Kep];
end
if length(vp)==1
	vp=[vp vp];
end
if length(Kt)==1
	Kt=[Kt Kt];
end

%----------------------------------------
%SIMULATE PERFUSION----------------------
%----------------------------------------

rt=exp(-Kep(1).*t);
rt_scar=exp(-Kep(2).*t);

if AIF_SIM
   %theoretical convolution
   %Concentration in tissue
   ct_RV=real(F1.*rt.*gammainc((t.*(1./b1-Kep(1))),a1)./(1-Kep(1).*b1).^a1);
   ct_LV=[zeros([1,delay_LV_RV]) ct_RV];
   %Concentration in scarr
   ct_scar=real(F1.*rt_scar.*gammainc((t.*(1./b1-Kep(2))),a1)./(1-Kep(2).*b1).^a1);
   ct_scar=[zeros([1,delay_LV_RV]) ct_scar];
end

%ct_RV=vp(1).*cat_RV +Cv_RV(1:length(cat_RV));
%ct_LV=vp(1).*cat_LV +Cv_LV(1:length(cat_LV));
%ct_scar=vp(2).*cat_LV +Cv_scar(1:length(cat_LV));

ct_RV=ct_RV(1:length(cat_RV));
ct_LV=ct_RV(1:length(cat_LV));
ct_scar=ct_scar(1:length(cat_LV));

%-----------------
ct_RV=[zeros([1,delay_perf]) ct_RV];
ct_LV=[zeros([1,delay_perf]) ct_LV];
ct_scar=[zeros([1,delay_perf]) ct_scar];


%----------------------------------------
%SIMULATE TISSUE WITH PERFUSION------------
%----------------------------------------
St=zeros([size(PH),Nt+1]);  %DISTORTED
PHt=zeros([size(PH),Nt+1]); %ORIGINAL VALUES
PHt(:,:,1)=PH;
St(:,:,1)=S0n;

%FOR ALL IMAGES IN PERFUSSION TIME
for ii=1:Nt
    if tmove>0
    	[Soi refM DefF]=motion_sim(S0,tmove,0);
	PHi=motion_sim(PH,tmove,1,refM,DefF);
	PHt(:,:,ii+1)=PHi;
        Bloodi=PHi==255;
        Blood2i=PHi==247;
	Blood3i=(PHi==128);
	LVi=(PHi==38);
	Scari=(PHi==26);

    else
	Soi=S0;
	Bloodi=Blood;
	Blood2i=Blood2;
	Blood3i=Blood3;
	LVi=LV;
	Scari=Scar;
	PHt(:,:,ii+1)=PH;
    end

T10i=-1./TR.*log(1-Soi./M0);

 Soi=Soi.*(1-(Bloodi+Blood2i+Blood3i+LVi+Scari))+...
    M0.*(1-exp(-TR.*(T10i+r1.*cat_LV(ii)))).*(Bloodi)+...
    M0.*(1-exp(-TR.*(T10i+r1.*cat_RV(ii)))).*(Blood2i)+...
    M0.*(1-exp(-TR.*(T10i+r1.*ct_LV(ii)))).*(LVi)+...
    M0.*(1-exp(-TR.*(T10i+r1.*ct_scar(ii)))).*(Scari)+...
    M0.*(1-exp(-TR.*(T10i+r1.*ct_LV(ii+1)))).*(Blood3i);  

St(:,:,ii+1)=distort_sg(Soi,sigma,blur,field, Knoise,ncorr,Ncoils);
if xmax>1
%St(:,:,ii+1)=min(St(:,:,ii+1),xmax);
end

BT(:,:,ii)=Bloodi;

end %END FOR


Iout=St;



%=======================================================
%========================FUNCIONS=======================
%========================================================
%========================================================
%PARSE INPUTS-------------------------------------------

function [Ncoil,ncorr,blur,field,parallel,tmove,xmax,AIF_SIM]  = parse_inputs(varargin)

Ncoil=1;
ncorr=1;
blur=1.5;
field=0;
tmove=0;
xmax=0;
parallel=0;
AIF_SIM=1;

flag=1;
for i = 1 : length(varargin)

   if flag==2
      flag=1;
     continue
    end

   if strcmp(varargin{i},'singlecoil')
      Ncoil = 1;
      flag=1;
   elseif strcmp(varargin{i},'multiplecoil')
      Ncoil=varargin{i+1};
      flag=1;
   elseif strcmp(varargin{i},'sense')
      parallel = 1;
      flag=1;
   elseif strcmp(varargin{i},'grappa')
      parallel = 2;
      flag=1;
   elseif strcmp(varargin{i},'No corr') 
      ncorr = 0;
      flag=1;
   elseif strcmp(varargin{i},'corr') 
      ncorr = 1;
      flag=1;
   elseif strcmp(varargin{i},'No PVE') 
      blur = 0;
      flag=1;
   elseif strcmp(varargin{i},'PVE') 
      blur=varargin{i+1};
      flag=2;
   elseif strcmp(varargin{i},'field')
      field=1;
      flag=1;
   elseif strcmp(varargin{i},'AIF_sim')
      AIF_SIM=1;
      flag=1;
   elseif strcmp(varargin{i},'AIF_file')
      AIF_SIM=0;
      flag=1;
   elseif strcmp(varargin{i},'tmove')
      tmove=varargin{i+1};
      flag=2;
   elseif  strcmp(varargin{i},'xmax')
      xmax=varargin{i+1};
      flag=2;
   end

end