 function pcorrected = stat_Zcorrected2D(stat, FWHM);
% pcorrected = stat_Zcorrected2D(stat, FWHM, img);
%
% stat_Zcorrected2D returns the corrected p-value of Z-statistic map with given 
% smoothing kernel sigma for 2D image 
%
% stat           : Given Z-statistics value
% FWHM      : smoothing Gaussian kernel size in full width at half maixmum (FWHM) in 
%                    terms of pixel size.
% pcorrected: corrected p-value
%
%
% References:
% [1] Worsley et al. 1995. A unified statistical approach for determining
% significant signals in location and scale sapce iamges of cerebral
% activation
% [2] Chung et al., 2012. Computational Neuroanatomy: The Methods, World Scientific Publishing
% 
%
% (C) 2012-2013  Moo K. Chung
% University of Wisconsin-Madison
% mkchung@wisc.edu
%
% Version 1. May 12, 2012 created
% Version 1.2 June 4, 2013 typo fixed
% Version 1.3 June 27, 2013 additional error fixed
   
n=size(stat);
t=reshape(stat, n(1)*n(2),1);
p=ones(n(1)*n(2),1);  % empty space is given p-value 1 by default

ind= find( abs(t)>1.64); % only Z-value above 1.64 or below -1.64 will be considered as a search area (SA)
area = length(ind); % area = total number of pixels in SA

tt=t(ind); % Z-value in SA
%rho2 = 4*log(2)/FWHM^2*tt.^2.*exp(-tt.^2/2)/(2*pi)^(3/2); % OLD CODE: code with typo in Version 1
tt=abs(tt); %Version 1.3
rho2 = 4*log(2)/FWHM^2*tt.*exp(-tt.^2/2)/(2*pi)^(3/2); % corrected code in Version 1.2

p(ind)= rho2*area;  % corrected p-value approximation. 

ind=find(p>0.1); %  if approximate pvalue > 0.1, let pavlue = 0.1 
p(ind)=0.1;
pcorrected = reshape(p, n(1), n(2));

subplot(2,2,4);; imagesc(pcorrected);  axis square; axis off; colorbar

maxc = max(max(max(pcorrected)), 0.1);  %max pvalue    
caxis([0 maxc]); 
title('p-value (corrected)')
