create_hrfpl_maps.m

function create_hrfpl_maps(Canonical, Derivative, Conditions, Mask)
%
% Creates HRF_PL images according to Henson et al., 2002.  Requires Matlab 
% 2017b or later for niftiread function. Otherwise change function for 
% reading and writing in volumes.
% Input: 
%       Canonical  -- file name of beta_????.img of canonical term
%       Derivative -- file name of beta_????.img of derivative term
%       Conditions -- cell with condition/session names as needed for output
%                       name, e.g.: cond  = {'PreStress', 'Stress', 'PostStress'}
%       Mask       -- file name of masking image. Optional (default is no
%                       mask)

% ----------------------------------------------------------------------
% 
% Authors: Immanuel Elbau and Philipp Sämann, MPI of Psychiatry- 2017
%
% Email: immanuel_elbau@psych.mpg.de or saemann@psych.mpg.de
% -------------------------------------------------------------------


%% Costants determined by optimizing fit of epmirical function with 'sigmoid_trans_optimization.m'

C(1,1:902629)   = -2;

D(1,1:902629)   = -0.7;    

Ycorr(1,1:902629)  =  0.25;

Xcorr(1,1:902629) = -0.5;

%% Load in beta images for canonical and derivative pictures

Info = niftiinfo(Canonical);

Canonical   = niftiread(Info);

Canonical = reshape(Canonical, [1 size(Canonical,1).*size(Canonical,2).*size(Canonical,3)]);

Info = niftiinfo(Derivative);

Derivative   = niftiread(Info);

Derivative = reshape(Derivative, [1 size(Derivative,1).*size(Derivative,2).*size(Derivative,3)]);


%% HRF-PL images according to Henson et al., 2002

% Calculate derivative / canonical ratio

HRFPL = Derivative./Canonical; 
        
% Perform sigmoid transformation

HRFPL  = 2 .* C ./(1+exp(D.*(HRFPL-Xcorr))) - C + Ycorr;
        
%% Apply amplitude mask.

% Load mask file and binarize it. If no mask file is entered as argument. 
% create an all ones mask (no masking)

if nargin < 4
    
    Mask = ones(1, size(Canonical,1).*size(Canonical,2).*size(Canonical,3));
    
else
    
    InfoMask = niftiinfo(Mask);

    Mask = niftiread(InfoMask);

    Mask = reshape(Mask, [1 size(Mask,1).*size(Mask,2).*size(Mask,3)]);
  
    % Binarize Mask
    
    Mask(Mask ~= 0 & ~isnan(Mask)) = 1;
    
    Mask(Mask~=1)     = 0;
    
    Mask(isnan(Mask)) = 0;
    
end

% Multiply HRF-PL image with mask

HRFPL = HRFPL .* Mask;

%% Write image

Output = ['HRFPL_map_',Conditions,'.nii'];

HRFPL = reshape(HRFPL, [Info.ImageSize]);

niftiwrite(HRFPL, Output, Info);

% The resulting images were smoothed as described in the SI Appendix (p.5) 
% with a 6 mm (FWHM) kernel
end