If you're just trying to clean things up then loops will help you get rid of a lot of repeated code when you are just changing a variable.
%read the asl (which is merged and motion corrected) to matlab
asl=load_nii('asl.nii.gz');
asl_img=double(asl.img);
%create two cell arrays one for the lable and the other for the control
l={asl_img(:,:,:,1); asl_img(:,:,:,3); asl_img(:,:,:,5); asl_img(:,:,:,7)};
c={asl_img(:,:,:,2); asl_img(:,:,:,4); asl_img(:,:,:,6); asl_img(:,:,:,8)};
%subtract both to create the subtracted asl images
dif_asl=cellfun(@minus,c,l, 'UniformOutput', false);
%split the cell aray thus each difrent-asl will be in a seprate cell
%{} to take the matrix out of the cell array, if we use (), example
%dif_asl1=dif_asl(1,1), it will remain as a cell aray and hence we will
%need the flowing command to create the mean comb1=cat(4,dif_asl1,dif_asl2);
dif_asl1=dif_asl{1,1}; % This section isn't necessary, just use indexing
dif_asl2=dif_asl{2,1};
dif_asl3=dif_asl{3,1};
dif_asl4=dif_asl{4,1};
%create the possipole combianations between the difference images
comb1=cat(4,dif_asl1,dif_asl2); % Run a loop: for k = 1:4; comb{k} = cat(4,dif_asl{k,1},dif_asl{k+1,1}) It's not perfect but it will simplify things
comb2=cat(4,dif_asl1,dif_asl3);
comb3=cat(4,dif_asl1,dif_asl4);
comb4=cat(4,dif_asl2,dif_asl3);
comb5=cat(4,dif_asl2,dif_asl4);
comb6=cat(4,dif_asl3,dif_asl4);
%mean of each combination
t2_comb1=nanmean(comb1,4); % Again, for loop: for k = 1:6; t2_comb{k} = nanmean(comb{k},4); end;
t2_comb2=nanmean(comb2,4);
t2_comb3=nanmean(comb3,4);
t2_comb4=nanmean(comb4,4);
t2_comb5=nanmean(comb5,4);
t2_comb6=nanmean(comb6,4);
%define the variable in order to calculate the CBF-map, because we use the values in ms, later on in the
%equation we multiply by 1000 to convert them to seconds and by 6000 to
%convert them from ml/g/s to ml/(100)g/min
y=0.9;
a=0.98; % beacause it is PASL
TI1=1000;
TI2=1800;
T1=1665;
%read the m0-map for calibration
S0=load_nii('m0s-m0map.nii');
S0=double(S0.img);
%create the CBF-maps from the created combinations
cbf1=asl;
cbf1.hdr.dime.dim([1,5])=[3,1]; % Perhaps the complexity of this structure is a default, but seems unnecessarily complicated since you're only defining two things
cbf1.img =(6000 * 1000 * y * exp(TI2/T1).*(t2_comb1))./(2*a*TI1 .*S0); % Again, for loop with indexing. Cell arrays can contain structures as well, cbf{k}
save_nii(cbf1,'cbf1.nii');
cbf2=asl;
cbf2.hdr.dime.dim([1,5])=[3,1];
cbf2.img =(6000 * 1000 * y * exp(TI2/T1).*(t2_comb2))./(2*a*TI1 .*S0);
save_nii(cbf2,'cbf2.nii');
cbf3=asl;
cbf3.hdr.dime.dim([1,5])=[3,1];
cbf3.img =(6000 * 1000 * y * exp(TI2/T1).*(t2_comb3))./(2*a*TI1 .*S0);
save_nii(cbf3,'cbf3.nii');
cbf4=asl;
cbf4.hdr.dime.dim([1,5])=[3,1];
cbf4.img =(6000 * 1000 * y * exp(TI2/T1).*(t2_comb4))./(2*a*TI1 .*S0);
save_nii(cbf4,'cbf4.nii');
cbf5=asl;
cbf5.hdr.dime.dim([1,5])=[3,1];
cbf5.img =(6000 * 1000 * y * exp(TI2/T1).*(t2_comb5))./(2*a*TI1 .*S0);
save_nii(cbf5,'cbf5.nii');
cbf6=asl;
cbf6.hdr.dime.dim([1,5])=[3,1];
cbf6.img =(6000 * 1000 * y * exp(TI2/T1).*(t2_comb6))./(2*a*TI1 .*S0);
save_nii(cbf6,'cbf6.nii');
%create the possible difference percentage between the map while keep the
%same dominator absolute(cbf1-cbf6)./((cbf1+cbf6)./2) where cbf1 is
%combination of pair 1 and 2 and cbf6 is combination 3 and 4
delta1.hdr=cbf1.hdr; % Last for loop here
delta1.img=abs(cbf1.img-cbf6.img)./abs((cbf1.img+cbf6.img)./2);
save_nii(delta1,'delta1.nii');
delta2.hdr=cbf1.hdr;
delta2.img=abs(cbf2.img-cbf5.img)./abs((cbf2.img+cbf5.img)./2);
save_nii(delta2,'delta2.nii');
delta3.hdr=cbf1.hdr;
delta3.img=abs(cbf3.img-cbf4.img)./abs((cbf3.img+cbf4.img)./2);
save_nii(delta3,'delta3.nii');
