function [outdir, prov] = spm_results_nidm(SPM,xSPM,TabDat)
% Export SPM stats results using the Neuroimaging Data Model (NIDM)
% FORMAT [outdir, prov] = spm_results_nidm(SPM,xSPM,TabDat)
% SPM     - structure containing analysis details (see spm_spm.m)
% xSPM    - structure containing inference details (see spm_getSPM.m)
% TabDat  - structure containing results details (see spm_list.m)
%
% outdir  - output directory
% prov    - provenance object (see spm_provenance.m)
%__________________________________________________________________________
% References:
% 
% Neuroimaging Data Model (NIDM):
%   http://nidm.nidash.org/
%
% PROV-DM: The PROV Data Model:
%   http://www.w3.org/TR/prov-dm/
%__________________________________________________________________________
% Copyright (C) 2013-2015 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_results_nidm.m 6468 2015-06-03 16:22:52Z guillaume $


%-Get input parameters, interactively if needed
%==========================================================================
if nargin < 1
    [SPM,xSPM] = spm_getSPM;
elseif nargin < 2
    if isstruct(SPM)
        xSPM = struct('swd',SPM.swd);
    else
        xSPM = struct('swd',spm_file('cpath',SPM));
    end
    [SPM,xSPM] = spm_getSPM(xSPM);
end
if nargin < 3
    % Consider Inf local maxima more than 0mm apart (i.e. all)
    TabDat = spm_list('Table',xSPM,Inf,0);
end

%-Options
%--------------------------------------------------------------------------
gz           = '.gz'; %-Compressed NIfTI {'.gz', ''}
coordsys     = nidm_conv('nidm:MNICoordinateSystem'); %-Assuming MNI space
NIDMversion  = '1.1.0';


%==========================================================================
%-Populate output directory
%==========================================================================
if ~exist(SPM.swd,'dir'), SPM.swd = pwd; end
outdir = fullfile(SPM.swd,'nidm');
outdir = spm_file(outdir,'uniquedir');
sts    = mkdir(outdir);
if ~sts, error('Cannot create directory "%s".',outdir); end

%-Design Matrix image (as png and csv)
%--------------------------------------------------------------------------
files.desimg = fullfile(outdir,'DesignMatrix.png');
DesMtx       = (SPM.xX.nKX + 1)*32;
ml           = floor(size(DesMtx,1)/size(DesMtx,2));
DesMtx       = reshape(repmat(DesMtx,ml,1),size(DesMtx,1),[]);
imwrite(DesMtx,gray(64),files.desimg,'png');
files.descsv = fullfile(outdir,'DesignMatrix.csv');
csvwrite(files.descsv,SPM.xX.xKXs.X);

%-Maximum Intensity Projection image (as png)
%--------------------------------------------------------------------------
files.mip    = fullfile(outdir,'MaximumIntensityProjection.png');
MIP          = spm_mip(xSPM.Z,xSPM.XYZmm,xSPM.M,xSPM.units);
imwrite(MIP,gray(64),files.mip,'png');

%-Beta images (as NIfTI)
%--------------------------------------------------------------------------
% for i=1:numel(SPM.Vbeta)
%     files.beta{i} = fullfile(outdir,[sprintf('Beta_%04d',i) '.nii' gz]);
%     img2nii(fullfile(xSPM.swd,SPM.Vbeta(i).fname), files.beta{i});
% end
files.beta = {};

%-SPM{.}, contrast, standard error and ESS images (as NIfTI)
%--------------------------------------------------------------------------
for i=1:numel(xSPM.Ic)
    if numel(xSPM.Ic) == 1, postfix = '';
    else                    postfix = sprintf('_%04d',i); end
    files.spm{i}  = fullfile(outdir,[xSPM.STAT 'Statistic' postfix '.nii' gz]);
    img2nii(fullfile(xSPM.swd,SPM.xCon(xSPM.Ic(i)).Vspm.fname), files.spm{i}, xSPM);
    if xSPM.STAT == 'T'
        files.con{i} = fullfile(outdir,['Contrast' postfix '.nii' gz]);
        img2nii(fullfile(xSPM.swd,SPM.xCon(xSPM.Ic(i)).Vcon.fname), files.con{i},...
            struct('STAT','con'));
        files.conse{i} = fullfile(outdir,['ContrastStandardError' postfix '.nii' gz]);
        Vc = SPM.xCon(xSPM.Ic(i)).c' * SPM.xX.Bcov * SPM.xCon(xSPM.Ic(i)).c;
    img2nii(fullfile(xSPM.swd,SPM.VResMS.fname), files.conse{i}, struct('fcn',@(x) sqrt(x*Vc)));
    elseif xSPM.STAT == 'F'
        files.ess{i} = fullfile(outdir,['ExtraSumOfSquares' postfix '.nii' gz]);
        img2nii(fullfile(xSPM.swd,SPM.xCon(xSPM.Ic(i)).Vcon.fname), files.ess{i});
    end
end

%-Thresholded SPM{.} image (as NIfTI)
%--------------------------------------------------------------------------
files.tspm = fullfile(outdir,['ExcursionSet.nii' gz]);
if ~isempty(gz), files.tspm = spm_file(files.tspm,'ext',''); end
spm_write_filtered(xSPM.Z,xSPM.XYZ,xSPM.DIM,xSPM.M,'',files.tspm);
if ~isempty(gz), gzip(files.tspm); spm_unlink(files.tspm); files.tspm = [files.tspm gz]; end

%-Residual Mean Squares image (as NIfTI)
%--------------------------------------------------------------------------
files.resms = fullfile(outdir,['ResidualMeanSquares.nii' gz]);
img2nii(fullfile(xSPM.swd,SPM.VResMS.fname), files.resms);

%-Resels per Voxel image (as NIfTI)
%--------------------------------------------------------------------------
files.rpv = fullfile(outdir,['ReselsPerVoxel.nii' gz]);
img2nii(fullfile(xSPM.swd,SPM.xVol.VRpv.fname), files.rpv);

%-Analysis mask image (as NIfTI)
%--------------------------------------------------------------------------
files.mask = fullfile(outdir,['Mask.nii' gz]);
img2nii(fullfile(xSPM.swd,SPM.VM.fname), files.mask);

%-Grand mean image (as NIfTI)
%--------------------------------------------------------------------------
files.grandmean = fullfile(outdir,'GrandMean.nii');
sf  = mean(SPM.xX.xKXs.X,1);
Vb  = SPM.Vbeta;
for i=1:numel(Vb), Vb(i).pinfo(1:2,:) = Vb(i).pinfo(1:2,:) * sf(i); end
Vgm = struct(...
    'fname',   files.grandmean,...
    'dim',     Vb(1).dim,...
    'dt',      [spm_type('float32') spm_platform('bigend')],...
    'mat',     Vb(1).mat,...
    'pinfo',   [1 0 0]',...
    'descrip', 'Grand Mean');
Vgm = spm_create_vol(Vgm);
Vgm.pinfo(1,1) = spm_add(Vb,Vgm);
Vgm = spm_create_vol(Vgm);
grandMeanMedian = spm_summarise(files.grandmean,SPM.VM.fname,@median);
if ~isempty(gz), gzip(files.grandmean); spm_unlink(files.grandmean); files.grandmean = [files.grandmean gz]; end

%-Explicit mask image (as NIfTI)
%--------------------------------------------------------------------------
if ~isempty(SPM.xM.VM)
    files.emask = fullfile(outdir,['CustomMask.nii' gz]);
    img2nii(fullfile(xSPM.swd,SPM.xM.VM.fname), files.emask);
else
    files.emask = '';
end

%-Clusters n-ary image (as NIfTI)
%--------------------------------------------------------------------------
files.clust = fullfile(outdir,['ClusterLabels.nii' gz]);
if ~isempty(gz), files.clust = spm_file(files.clust,'ext',''); end
Z   = spm_clusters(xSPM.XYZ);
idx = find(~cellfun(@isempty,{TabDat.dat{:,5}}));
n   = zeros(1,numel(idx));
for i=1:numel(idx)
    [unused,j] = spm_XYZreg('NearestXYZ',TabDat.dat{idx(i),12}',xSPM.XYZmm);
    n(i) = Z(j);
end
n(n) = 1:numel(n);
if max(Z) ~= numel(idx)
    warning('Small Volume Correction not handled yet.');
    n(numel(idx)+1:max(Z)) = 0;
end
Z    = n(Z);
spm_write_filtered(Z,xSPM.XYZ,xSPM.DIM,xSPM.M,'',files.clust);
if ~isempty(gz), gzip(files.clust); spm_unlink(files.clust); files.clust = [files.clust gz]; end

%-Display mask images (as NIfTI)
%--------------------------------------------------------------------------
for i=1:numel(xSPM.Im)
    files.dmask{i} = fullfile(outdir,[sprintf('DisplayMask_%04d.nii',i) gz]);
    if isnumeric(xSPM.Im)
        um = spm_u(xSPM.pm,[SPM.xCon(xSPM.Im(i)).eidf,SPM.xX.erdf],...
            SPM.xCon(xSPM.Im(i)).STAT);
        if ~xSPM.Ex, fcn = @(x) x > um;
        else         fcn = @(x) x <= um; end
        img2nii(SPM.xCon(xSPM.Im(i)).Vspm.fname, files.dmask{i}, struct('fcn',fcn));
    else
        if ~xSPM.Ex, fcn = @(x) x~=0 & ~isnan(x);
        else         fcn = @(x) ~(x~=0 & ~isnan(x)); end
        img2nii(xSPM.Im{i}, files.dmask{i}, struct('fcn',fcn));
    end
end
if numel(xSPM.Im) == 0, files.dmask = {}; end

%-SVC Mask (as NIfTI)
%--------------------------------------------------------------------------
if strcmp(TabDat.tit,'p-values adjusted for search volume')
    files.svcmask = '';
elseif strncmp(TabDat.tit,'search volume: ',15)
    warning('Small Volume Correction not handled yet.'); % see spm_VOI.m
    % '%0.1fmm sphere at [%.0f,%.0f,%.0f]'
    % '%0.1f x %0.1f x %0.1f mm box at [%.0f,%.0f,%.0f]'
    % 'image mask: %s'
    files.svcmask = '';
else
    warning('Unable to retrieve search volume details: assuming whole brain search.');
    files.svcmask = '';
end

%-Search Space mask image (as NIfTI)
%--------------------------------------------------------------------------
files.searchspace = fullfile(outdir,['SearchSpaceMask.nii' gz]);
img2nii(fullfile(xSPM.swd,SPM.VM.fname), files.searchspace);


%==========================================================================
%-                          D A T A   M O D E L
%==========================================================================

clear coordspace originalfile isHumanReadable

niifmt = {'image/nifti','xsd:string'};
isHumanReadable(true);

pp = spm_provenance;

%-Namespaces
%--------------------------------------------------------------------------
pp.add_namespace('nidm','http://purl.org/nidash/nidm#');
pp.add_namespace('niiri','http://iri.nidash.org/');
pp.add_namespace('spm','http://purl.org/nidash/spm#');
pp.add_namespace('neurolex','http://neurolex.org/wiki/');
pp.add_namespace('crypto','http://id.loc.gov/vocabulary/preservation/cryptographicHashFunctions#');
pp.add_namespace('dct','http://purl.org/dc/terms/');
pp.add_namespace('nfo','http://www.semanticdesktop.org/ontologies/2007/03/22/nfo#');
pp.add_namespace('dc','http://purl.org/dc/elements/1.1/');
pp.add_namespace('dctype','http://purl.org/dc/dcmitype/');
pp.add_namespace('obo','http://purl.obolibrary.org/obo/');

%-Provenance
%--------------------------------------------------------------------------
idResults = getid('niiri:spm_results_id',isHumanReadable);
pp.entity(idResults,{...
  'prov:type','prov:Bundle',...
  'prov:type',nidm_conv('nidm:Results'),...
  'prov:label','NIDM-Results',...
  nidm_conv('nidm:version'),{NIDMversion,'xsd:string'},...
  });
pp.wasGeneratedBy(idResults,'-',now);

p = spm_provenance;

%-Agent: SPM
%--------------------------------------------------------------------------
[V,R] = spm('Ver');
idSoftware = getid('niiri:software_id',isHumanReadable);
p.agent(idSoftware,{...
    'prov:type',nidm_conv('neurolex:SPM'),...
    'prov:type','prov:SoftwareAgent',...
    'prov:label',{'SPM','xsd:string'},...
    nidm_conv('nidm:softwareVersion'),{[V(4:end) '.' R],'xsd:string'},...
    });

%-Entity: Coordinate Space
%--------------------------------------------------------------------------
id_data_coordspace = coordspace(p,xSPM.M,xSPM.DIM,xSPM.units,coordsys,1);

%-Entity: Image Data
%--------------------------------------------------------------------------
if isfield(SPM,'Sess')
    extra_fields = {...
        nidm_conv('nidm:grandMeanScaling'),{'true','xsd:boolean'},...
        nidm_conv('nidm:targetIntensity'),{SPM.xGX.GM,'xsd:float'},...
        };
else
    extra_fields = {...
        nidm_conv('nidm:grandMeanScaling'),{'false','xsd:boolean'},...
        };
end
idData = getid('niiri:data_id',isHumanReadable);
p.entity(idData,{...
    'prov:type','prov:Collection',...
    'prov:type',nidm_conv('nidm:Data'),...
    'prov:label',{'Data','xsd:string'},...
    extra_fields{:}});

%-Entity: Drift Model
%--------------------------------------------------------------------------
if isfield(SPM,'Sess') && isfield(SPM.xX,'K')
    idDriftModel = getid('niiri:drift_model_id',isHumanReadable);
    
    p.entity(idDriftModel,{...
        'prov:type',nidm_conv('spm:DCTDriftModel'),...
        'prov:label','SPM''s DCT Drift Model',...
        nidm_conv('spm:driftCutoffPeriod'),{SPM.xX.K(1).HParam,'xsd:float'},...
        });
    extra_fields_drift = {nidm_conv('nidm:hasDriftModel'),idDriftModel};
else
    extra_fields_drift = {};
end

%-Entity: Design Matrix
%--------------------------------------------------------------------------
idDesignMatrix = getid('niiri:design_matrix_id',isHumanReadable);
idDesignMatrixImage = getid('niiri:design_matrix_png_id',isHumanReadable);

extra_fields_basis_set = {};
if isfield(SPM,'xBF')
    switch SPM.xBF.name
        case 'hrf'
            extra_fields_basis_set = ...
                {nidm_conv('nidm:hasHRFBasis'),nidm_conv('spm:GammaDifferenceHRF')};
        case 'hrf (with time derivative)'
            extra_fields_basis_set = {...
                nidm_conv('nidm:hasHRFBasis'),nidm_conv('spm:GammaDifferenceHRF'),...
                nidm_conv('nidm:hasHRFBasis'),nidm_conv('spm:TemporalDerivative')};
        case 'hrf (with time and dispersion derivatives)'
            extra_fields_basis_set = {...
                nidm_conv('nidm:hasHRFBasis'),nidm_conv('spm:GammaDifferenceHRF'),...
                nidm_conv('nidm:hasHRFBasis'),nidm_conv('spm:TemporalDerivative'),...
                nidm_conv('nidm:hasHRFBasis'),nidm_conv('spm:DispersionDerivative')};
        case {'Fourier set','Fourier set (Hanning)','Gamma functions','Finite Impulse Response'}
            warning('Not implemented yet.');
        otherwise
            warning('Unknown basis set.');
    end
end

p.entity(idDesignMatrix,{...
    'prov:type',nidm_conv('nidm:DesignMatrix'),...
    'prov:location',{uri(spm_file(files.descsv,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.descsv,'filename'),'xsd:string'},...
    'dct:format',{'text/csv','xsd:string'},...
    'dc:description',idDesignMatrixImage,...
    'prov:label',{'Design Matrix','xsd:string'},...
    nidm_conv('nidm:regressorNames'),{SPM.xX.name,'xsd:string'},...
    extra_fields_drift{:},...
    extra_fields_basis_set{:}
    });
    %nidm_conv('nidm:hasfMRIDesign'),nidm_conv('nidm:BlockDesign'),...

p.entity(idDesignMatrixImage,{...
    'prov:type','dctype:Image',...
    'prov:location',{uri(spm_file(files.desimg,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.desimg,'filename'),'xsd:string'},...
    'dct:format',{'image/png','xsd:string'},...
    });

%-Entity: Explicit Mask
%--------------------------------------------------------------------------
if ~isempty(SPM.xM.VM)
    if ~spm_check_orientations(struct('dim',{xSPM.DIM',SPM.xM.VM.dim},...
            'mat',{xSPM.M,SPM.xM.VM.mat}),false)
        id_emask_coordspace = coordspace(p,SPM.xM.VM.mat,SPM.xM.VM.dim',...
            xSPM.units,coordsys);
    end
    idMask2 = getid('niiri:mask_id_2',isHumanReadable);
    p.entity(idMask2,{...
        'prov:type','nidm:CustomMaskMap',...
        'prov:location',{uri(spm_file(files.emask,'cpath')),'xsd:anyURI'},...
        'nfo:fileName',{spm_file(files.emask,'filename'),'xsd:string'},...
        'dct:format',niifmt,...
        'prov:label',{'Custom Mask','xsd:string'},...
        nidm_conv('nidm:inCoordinateSpace'),id_emask_coordspace,...
        'crypto:sha512',{sha512sum(spm_file(files.emask,'cpath')),'xsd:string'},...
        });
    id = originalfile(p,files.emask,idMask2,'nidm:CustomMaskMap');
    p.wasDerivedFrom(idMask2,id);
end

%-Entity: Error Model
%--------------------------------------------------------------------------
if isfield(SPM.xVi,'form')
    if strcmp(SPM.xVi.form,'i.i.d')
        extra_fields_NM = {...
            nidm_conv('nidm:hasErrorDependence'),nidm_conv('nidm:IndependentError'),...
            nidm_conv('nidm:errorVarianceHomogeneous'),{'true','xsd:boolean'},...
            };
        extra_fields_PE = {
            nidm_conv('nidm:withEstimationMethod'),nidm_conv('stato:OLS'),...
            };
    else
        extra_fields_NM = {...
            nidm_conv('nidm:hasErrorDependence'),nidm_conv('stato:ToeplitzCovarianceStructure'),...
            nidm_conv('nidm:dependenceSpatialModel'),nidm_conv('nidm:SpatiallyGlobalModel'),...
            nidm_conv('nidm:errorVarianceHomogeneous'),{'true','xsd:boolean'},...
            nidm_conv('nidm:varianceSpatialModel'),nidm_conv('nidm:SpatiallyLocalModel'),...
            };
        extra_fields_PE = {
            nidm_conv('nidm:withEstimationMethod'),nidm_conv('stato:GLS'),...
            };
    end
else
    if ~isfield(SPM.xVi,'Vi') || numel(SPM.xVi.Vi) == 1 % assume it's identity
        extra_fields_NM = {...
            nidm_conv('nidm:hasErrorDependence'),nidm_conv('nidm:IndependentError'),...
            nidm_conv('nidm:errorVarianceHomogeneous'),{'true','xsd:boolean'},...
            };
        extra_fields_PE = {
            nidm_conv('nidm:withEstimationMethod'),nidm_conv('stato:OLS'),...
            };
    else
        extra_fields_NM = {...
            nidm_conv('nidm:hasErrorDependence'),'nidm:ArbitrarilyCorrelatedError',...
            'nidm:dependenceSpatialModel','nidm:SpatiallyGlobalModel',...
            'nidm:errorVarianceHomogeneous',{'false','xsd:boolean'},...
            'nidm:varianceSpatialModel','nidm:SpatiallyLocalModel',...
            };
        extra_fields_PE = {
            nidm_conv('nidm:withEstimationMethod'),nidm_conv('stato:GLS'),...
            };
    end
end
idErrorModel = getid('niiri:error_model_id',isHumanReadable);
p.entity(idErrorModel,{...
    'prov:type',nidm_conv('nidm:ErrorModel'),...
    nidm_conv('nidm:hasErrorDistribution'),nidm_conv('nidm:GaussianDistribution'),...
    extra_fields_NM{:}});

%-Activity: Model Parameters Estimation
%==========================================================================
idModelPE = getid('niiri:model_pe_id',isHumanReadable);
p.activity(idModelPE,{...
    'prov:type',nidm_conv('nidm:ModelParametersEstimation'),...
    'prov:label','Model parameters estimation',...
    extra_fields_PE{:}});
p.wasAssociatedWith(idModelPE, idSoftware);
p.used(idModelPE, idDesignMatrix);
p.used(idModelPE, idData);
p.used(idModelPE, idErrorModel);
if ~isempty(SPM.xM.VM)
    p.used(idModelPE, idMask2);
end

%-Entity: Mask Map
%--------------------------------------------------------------------------
idMask1 = getid('niiri:mask_id_1',isHumanReadable);
p.entity(idMask1,{...
    'prov:type',nidm_conv('nidm:MaskMap'),...
    'prov:location',{uri(spm_file(files.mask,'cpath')),'xsd:anyURI'},...
    nidm_conv('nidm:isUserDefined'),{'false','xsd:boolean'},...
    'nfo:fileName',{spm_file(files.mask,'filename'),'xsd:string'},...
    'dct:format',niifmt,...
    'prov:label',{'Mask','xsd:string'},...
    nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
    'crypto:sha512',{sha512sum(spm_file(files.mask,'cpath')),'xsd:string'},...
    });
id = originalfile(p,SPM.VM.fname,idMask1,nidm_conv('nidm:MaskMap'));
p.wasDerivedFrom(idMask1,id);
p.wasGeneratedBy(idMask1, idModelPE);

%-Entity: Grand Mean Map
%--------------------------------------------------------------------------
idGrandMean = getid('niiri:grand_mean_map_id',isHumanReadable);
p.entity(idGrandMean,{...
    'prov:type',nidm_conv('nidm:GrandMeanMap'),...
    'prov:location',{uri(spm_file(files.grandmean,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.grandmean,'filename'),'xsd:string'},...
    'dct:format',niifmt,...
    'prov:label',{'Grand Mean Map','xsd:string'},...
    nidm_conv('nidm:maskedMedian'),{grandMeanMedian,'xsd:float'},...
    nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
    'crypto:sha512',{sha512sum(spm_file(files.grandmean,'cpath')),'xsd:string'},...
    });
p.wasGeneratedBy(idGrandMean, idModelPE);

%-Entity: Beta Maps
%--------------------------------------------------------------------------
idBeta = cell(1,numel(SPM.Vbeta));
for i=1:numel(SPM.Vbeta)
    if ~isempty(files.beta)
        extra_fields = {...
            'prov:location',{uri(files.beta{i}),'xsd:anyURI'},...
            'nfo:fileName',{spm_file(files.beta{i},'filename'),'xsd:string'},...
            'dct:format',niifmt,...
            'crypto:sha512',{sha512sum(files.beta{i}),'xsd:string'},...
        };
    else
        extra_fields = {};
    end
    idBeta{i} = getid(sprintf('niiri:beta_map_id_%d',i),isHumanReadable);
    p.entity(idBeta{i},{...
        'prov:type',nidm_conv('nidm:ParameterEstimateMap'),...
        'prov:label',{sprintf('Beta Map %d',i),'xsd:string'},...
        nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
        extra_fields{:},...
    });
    id = originalfile(p,SPM.Vbeta(i).fname,idBeta{i},nidm_conv('nidm:ParameterEstimateMap'));
    p.wasDerivedFrom(idBeta{i},id);
    p.wasGeneratedBy(idBeta{i}, idModelPE);
end

%-Entity: ResMS Map
%--------------------------------------------------------------------------
idResMS = getid('niiri:residual_mean_squares_map_id',isHumanReadable);
p.entity(idResMS,{...
    'prov:type',nidm_conv('nidm:ResidualMeanSquaresMap'),...
    'prov:location',{uri(spm_file(files.resms,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.resms,'filename'),'xsd:string'},...
    'dct:format',niifmt,...
    'prov:label',{'Residual Mean Squares Map','xsd:string'},...
    nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
    'crypto:sha512',{sha512sum(spm_file(files.resms,'cpath')),'xsd:string'},...
    });
id = originalfile(p,SPM.VResMS.fname,idResMS,nidm_conv('nidm:ResidualMeanSquaresMap'));
p.wasDerivedFrom(idResMS,id);
p.wasGeneratedBy(idResMS, idModelPE);

%-Entity: RPV Map
%--------------------------------------------------------------------------
idRPV = getid('niiri:resels_per_voxel_map_id',isHumanReadable);
p.entity(idRPV,{...
    'prov:type',nidm_conv('nidm:ReselsPerVoxelMap'),...
    'prov:location',{uri(spm_file(files.rpv,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.rpv,'filename'),'xsd:string'},...
    'dct:format',niifmt,...
    'prov:label',{'Resels per Voxel Map','xsd:string'},...
    nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
    'crypto:sha512',{sha512sum(spm_file(files.rpv,'cpath')),'xsd:string'},...
    });
id = originalfile(p,SPM.xVol.VRpv.fname,idRPV,nidm_conv('nidm:ReselsPerVoxelMap'));
p.wasDerivedFrom(idRPV,id);
p.wasGeneratedBy(idRPV, idModelPE);

%-Activity: Contrast Estimation
%==========================================================================
for c=1:numel(xSPM.Ic)
    if numel(xSPM.Ic) == 1, postfix = '';
    else                    postfix = sprintf('_%d',c); end
    
    idConVec = getid(['niiri:contrast_id' postfix],isHumanReadable);
    p.entity(idConVec,{...
        'prov:type',nidm_conv('stato:ContrastWeightMatrix'),...
        nidm_conv('nidm:statisticType'),nidm_conv(['stato:' xSPM.STAT 'Statistic']),...
        nidm_conv('nidm:contrastName'),{SPM.xCon(xSPM.Ic(c)).name,'xsd:string'},... %esc
        'prov:label',{['Contrast: ' SPM.xCon(xSPM.Ic(c)).name],'xsd:string'},... %esc
        'prov:value',{SPM.xCon(xSPM.Ic(c)).c','xsd:string'},...
        });

    idConEst = getid(['niiri:contrast_estimation_id' postfix],isHumanReadable);
    p.activity(idConEst,{...
        'prov:type',nidm_conv('nidm:ContrastEstimation'),...
        'prov:label',['Contrast estimation' strrep(postfix,'_',' ')],...
        });
    p.wasAssociatedWith(idConEst, idSoftware);
    p.used(idConEst, idMask1);
    p.used(idConEst, idResMS);
    p.used(idConEst, idDesignMatrix);
    p.used(idConEst,idConVec);
    for i=1:numel(SPM.Vbeta)
        p.used(idConEst, idBeta{i});
    end
    
    idSPM{c} = getid(['niiri:statistic_map_id' postfix],isHumanReadable);
    p.entity(idSPM{c},{...
        'prov:type',nidm_conv('nidm:StatisticMap'),...
        'prov:location',{uri(spm_file(files.spm{c},'cpath')),'xsd:anyURI'},...
        'nfo:fileName',{spm_file(files.spm{c},'filename'),'xsd:string'},...
        'dct:format',niifmt,...
        'prov:label',{['Statistic Map: ' SPM.xCon(xSPM.Ic(c)).name],'xsd:string'},... %esc
        nidm_conv('nidm:statisticType'),nidm_conv(['stato:' xSPM.STAT 'Statistic']),...
        nidm_conv('nidm:contrastName'),{SPM.xCon(xSPM.Ic(c)).name,'xsd:string'},... %esc
        nidm_conv('nidm:errorDegreesOfFreedom'),{xSPM.df(2),'xsd:float'},...
        nidm_conv('nidm:effectDegreesOfFreedom'),{xSPM.df(1),'xsd:float'},...
        nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
        'crypto:sha512',{sha512sum(files.spm{c}),'xsd:string'},...
        });
    id = originalfile(p,SPM.xCon(xSPM.Ic(c)).Vspm.fname,idSPM{c},nidm_conv('nidm:StatisticMap'));
    p.wasDerivedFrom(idSPM{c},id);
    p.wasGeneratedBy(idSPM{c},idConEst);
    
    if xSPM.STAT == 'T'
        idContrast = getid(['niiri:contrast_map_id' postfix],isHumanReadable);
        p.entity(idContrast,{...
            'prov:type',nidm_conv('nidm:ContrastMap'),...
            'prov:location',{uri(spm_file(files.con{c},'cpath')),'xsd:anyURI'},...
            'nfo:fileName',{spm_file(files.con{c},'filename'),'xsd:string'},...
            'dct:format',niifmt,...
            'prov:label',{['Contrast Map: ' SPM.xCon(xSPM.Ic(c)).name],'xsd:string'},... %esc
            nidm_conv('nidm:contrastName'),{SPM.xCon(xSPM.Ic(c)).name,'xsd:string'},... %esc
            nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
            'crypto:sha512',{sha512sum(spm_file(files.con{c},'cpath')),'xsd:string'},...
            });
        id = originalfile(p,SPM.xCon(xSPM.Ic(c)).Vcon.fname,idContrast,nidm_conv('nidm:ContrastMap'));
        p.wasDerivedFrom(idContrast,id);
        p.wasGeneratedBy(idContrast,idConEst);
        
        idSE = getid(['niiri:contrast_standard_error_map_id' postfix],isHumanReadable);
        p.entity(idSE,{...
            'prov:type',nidm_conv('nidm:ContrastStandardErrorMap'),...
            'prov:location',{uri(spm_file(files.conse{c},'cpath')),'xsd:anyURI'},...
            'nfo:fileName',{spm_file(files.conse{c},'filename'),'xsd:string'},...
            'dct:format',niifmt,...
            'prov:label',{'Contrast Standard Error Map','xsd:string'},...
            nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
            'crypto:sha512',{sha512sum(spm_file(files.conse{c},'cpath')),'xsd:string'},...
            });
        p.wasGeneratedBy(idSE,idConEst);
    end
    if xSPM.STAT == 'F'
        idESS = getid(['niiri:contrast_extra_sum_of_squares_id' postfix],isHumanReadable);
        p.entity(idESS,{...
            'prov:type','nidm:ContrastExtraSumOfSquaresMap',...
            'prov:location',{uri(spm_file(files.ess{c},'cpath')),'xsd:anyURI'},...
            'nfo:fileName',{spm_file(files.ess{c},'filename'),'xsd:string'},...
            'dct:format',niifmt,...
            'prov:label',{'Contrast Extra Sum of Squares Map','xsd:string'},...
            nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
            'crypto:sha512',{sha512sum(spm_file(files.ess{c},'cpath')),'xsd:string'},...
            });
        id = originalfile(p,SPM.xCon(xSPM.Ic(c)).Vcon.fname,idESS,'nidm:ContrastExtraSumOfSquaresMap');
        p.wasDerivedFrom(idESS,id);
        p.wasGeneratedBy(idESS,idConEst);
    end
end

%-Entity: Height Threshold
%--------------------------------------------------------------------------
thresh(1).type  = nidm_conv('obo:Statistic');
thresh(1).label = 'Height Threshold';
thresh(1).value = xSPM.u; % TabDat.ftr{1,2}(1)
thresh(2).type  = nidm_conv('nidm:pValueUncorrectedClass');
thresh(2).label = 'Height Threshold';
thresh(2).value = TabDat.ftr{1,2}(2);
thresh(3).type  = nidm_conv('obo:pValueFWER');
thresh(3).label = 'Height Threshold';
thresh(3).value = TabDat.ftr{1,2}(3);
td = regexp(xSPM.thresDesc,'p\D?(?<u>[\.\d]+) \((?<thresDesc>\S+)\)','names');
if isempty(td)
    td = regexp(xSPM.thresDesc,'\w=(?<u>[\.\d]+)','names');
    if ~isempty(td)
        thresh_order = 1:3; % Statistic
        thresh_desc  = sprintf(': %s=%f)',xSPM.STAT,xSPM.u);
    else
        warning('Unkwnown threshold type.');
        thresh_order = 1:3; % unknown
        thresh_desc  = '';
    end
else
    switch td.thresDesc
        case 'FWE'
            thresh_order = [3 1 2]; % FWE
            thresh_desc  = sprintf(': p<%f (FWE)',TabDat.ftr{1,2}(3));
        case 'unc.'
            thresh_order = [2 1 3]; % uncorrected
            thresh_desc  = sprintf(': p<%f (unc.)',TabDat.ftr{1,2}(2));
        case 'FDR'
            warning('FDR not handled.');
            thresh_order = 1:3; % unknown
            thresh_desc  = '';
        otherwise
            warning('Unkwnown threshold type.');
            thresh_order = 1:3; % unknown
            thresh_desc  = '';
    end
end
thresh = thresh(thresh_order);
thresh(1).label = [thresh(1).label thresh_desc];
idHeightThresh = getid('niiri:height_threshold_id',isHumanReadable);
idHeightThresh2 = getid('niiri:height_threshold_id_2',isHumanReadable);
idHeightThresh3 = getid('niiri:height_threshold_id_3',isHumanReadable);
p.entity(idHeightThresh,{...
    'prov:type',nidm_conv('nidm:HeightThreshold'),...
    'prov:label',{thresh(1).label,'xsd:string'},... %esc
    'prov:value',{thresh(1).value,'xsd:float'},...
    nidm_conv('nidm:equivalentThreshold'),idHeightThresh2,...
    nidm_conv('nidm:equivalentThreshold'),idHeightThresh3,...
    });
p.entity(idHeightThresh2,{...
    'prov:type',nidm_conv('nidm:HeightThreshold'),...
    'prov:type',thresh(2).type,...
    'prov:label',{thresh(2).label,'xsd:string'},... %esc
    'prov:value',{thresh(2).value,'xsd:float'},...
    });

p.entity(idHeightThresh3,{...
    'prov:type',nidm_conv('nidm:HeightThreshold'),...
    'prov:type',thresh(3).type,...
    'prov:label',{thresh(3).label,'xsd:string'},... %esc
    'prov:value',{thresh(3).value,'xsd:float'},...
    });

%-Entity: Extent Threshold
%--------------------------------------------------------------------------
if spm_get_defaults('stats.rft.nonstat')
    warning('Non-stationary RFT results not handled yet.');
end
V2R = 1 / prod(xSPM.FWHM);

if TabDat.ftr{2,2}(1) > 0
    kk = [TabDat.ftr{2,2}(2) TabDat.ftr{2,2}(3)];
else
    kk = [1 1];
end
idExtentThresh = getid('niiri:extent_threshold_id',isHumanReadable);
idExtentThresh2 = getid('niiri:extent_threshold_id_2',isHumanReadable);
idExtentThresh3 = getid('niiri:extent_threshold_id_3',isHumanReadable);
p.entity(idExtentThresh,{...
    'prov:type',nidm_conv('nidm:ExtentThreshold'),...
    'prov:type',nidm_conv('obo:Statistic'),...
    'prov:label',{['Extent Threshold: k=' num2str(TabDat.ftr{2,2}(1))],'xsd:string'},...
    nidm_conv('nidm:clusterSizeInVoxels'),{TabDat.ftr{2,2}(1),'xsd:int'},... % xSPM.k
    nidm_conv('nidm:clusterSizeInResels'),{TabDat.ftr{2,2}(1)*V2R,'xsd:float'},...
    nidm_conv('nidm:equivalentThreshold'),idExtentThresh2,...
    nidm_conv('nidm:equivalentThreshold'),idExtentThresh3,...
    });
p.entity(idExtentThresh2,{...
    'prov:type',nidm_conv('nidm:ExtentThreshold'),...
    'prov:type',nidm_conv('obo:pValueFWER'),...
    'prov:label',{'Extent Threshold','xsd:string'},...
    'prov:value',{kk(2),'xsd:float'},...
    });
p.entity(idExtentThresh3,{...
    'prov:type',nidm_conv('nidm:ExtentThreshold'),...
    'prov:type',nidm_conv('nidm:pValueUncorrectedClass'),...
    'prov:label',{'Extent Threshold','xsd:string'},...
    'prov:value',{kk(1),'xsd:float'},...
    });

%-Entity: Peak Definition Criteria
%--------------------------------------------------------------------------
% TabDat.str = 'table shows %d local maxima more than %.1fmm apart'
maxNumberOfPeaksPerCluster = spm_get_defaults('stats.results.volume.nbmax');
minDistanceBetweenPeaks = spm_get_defaults('stats.results.volume.distmin');
% clusterConnectivityCriterion = 18;
idPeakDefCrit = getid('niiri:peak_definition_criteria_id',isHumanReadable);
p.entity(idPeakDefCrit,{...
    'prov:type',nidm_conv('nidm:PeakDefinitionCriteria'),...
    'prov:label',{'Peak Definition Criteria','xsd:string'},...
    nidm_conv('nidm:maxNumberOfPeaksPerCluster'),{maxNumberOfPeaksPerCluster,'xsd:int'},...
    nidm_conv('nidm:minDistanceBetweenPeaks'),{minDistanceBetweenPeaks,'xsd:float'},...
    });

%-Entity: Cluster Definition Criteria
%--------------------------------------------------------------------------
clusterConnectivityCriterion = 18; % see spm_max.m
idClusterDefCrit = getid('niiri:cluster_definition_criteria_id',isHumanReadable);
p.entity(idClusterDefCrit,{...
    'prov:type',nidm_conv('nidm:ClusterDefinitionCriteria'),...
    'prov:label',{sprintf('Cluster Connectivity Criterion: %d',clusterConnectivityCriterion),'xsd:string'},...
    nidm_conv('nidm:hasConnectivityCriterion'),nidm_conv(sprintf('nidm:voxel%dconnected',clusterConnectivityCriterion)),...
    });


%-Activity: Inference
%==========================================================================
if numel(xSPM.Ic) == 1
    st = {'prov:type',nidm_conv('nidm:Inference'), ...
          nidm_conv('nidm:hasAlternativeHypothesis'),nidm_conv('nidm:OneTailedTest'),...
          'prov:label','Inference'};
else
    if xSPM.n == 1
        st = {'prov:type',nidm_conv('nidm:ConjunctionInference'), ...
              'prov:label','Conjunction Inference'};
    else
        st = {'prov:type','nidm:kConjunctionInference', ...
              'prov:label','k-Conjunction Inference', ...
              nidm_conv('spm:partialConjunctionDegree'),xSPM.n};
    end
end
idInference = getid('niiri:inference_id',isHumanReadable);
p.activity(idInference,st);
p.wasAssociatedWith(idInference, idSoftware);
p.used(idInference, idHeightThresh);
p.used(idInference, idExtentThresh);
for c=1:numel(xSPM.Ic)
    if numel(xSPM.Ic) == 1, postfix = '';
    else                    postfix = sprintf('_%d',c); end
    p.used(idInference, idSPM{c});
end
p.used(idInference, idRPV);
p.used(idInference, idMask1);
p.used(idInference, idPeakDefCrit);
p.used(idInference, idClusterDefCrit);

%-Entity: Display Mask Maps
%--------------------------------------------------------------------------
for i=1:numel(files.dmask)
    V = spm_vol(files.dmask{i});
    if ~spm_check_orientations(struct('dim',{xSPM.DIM',V.dim},...
            'mat',{xSPM.M,V.mat}),false)
        currCoordSpace = coordspace(p,V.mat,V.dim',xSPM.units,coordsys);
    else
        currCoordSpace = id_data_coordspace;
    end
    
    if numel(files.dmask) == 1, postfix = '';
    else                    postfix = sprintf('_%d',i); end
    idDMask = getid(['niiri:display_mask_map_id' postfix],isHumanReadable);
    p.entity(idDMask,{...
        'prov:type',nidm_conv('nidm:DisplayMaskMap'),...
        'prov:location',{uri(spm_file(files.dmask{i},'cpath')),'xsd:anyURI'},...
        'nfo:fileName',{spm_file(files.dmask{i},'filename'),'xsd:string'},...
        'dct:format',niifmt,...
        'prov:label',{'Display Mask Map','xsd:string'},...
        nidm_conv('nidm:inCoordinateSpace'),currCoordSpace,...
        'crypto:sha512',{sha512sum(spm_file(files.dmask{i},'cpath')),'xsd:string'},...
        });
    id = originalfile(p,files.dmask{i},idDMask,nidm_conv('nidm:DisplayMaskMap'));
    p.wasDerivedFrom(idDMask,id);
    p.used(idInference, idDMask);
end

%-Entity: SVC Mask Map
%--------------------------------------------------------------------------
if ~isempty(files.svcmask)
    V = spm_vol(files.svcmask);
    if ~spm_check_orientations(struct('dim',{xSPM.DIM',V.dim},...
            'mat',{xSPM.M,V.mat}),false)
        currCoordSpace = coordspace(p,V.mat,V.dim',xSPM.units,coordsys);
    else
        currCoordSpace = id_data_coordspace;
    end
    idSVC = getid('niiri:sub_volume_id',isHumanReadable);
    p.entity(idSVC,{...
        'prov:type','nidm:SubVolumeMap',...
        'prov:location',{uri(spm_file(files.svcmask,'cpath')),'xsd:anyURI'},...
        'nfo:fileName',{spm_file(files.svcmask,'filename'),'xsd:string'},...
        'prov:label',{'Sub-volume','xsd:string'},...
        nidm_conv('nidm:inCoordinateSpace'),currCoordSpace,...
        'crypto:sha512',{sha512sum(spm_file(files.svcmask,'cpath')),'xsd:string'},...
        });
    id = originalfile(p,files.svcmask,idSVC,'nidm:SubVolumeMap');
    p.wasDerivedFrom(idSVC,id);
    p.used(idInference, idSVC);
end

%-Entity: Search Space
%--------------------------------------------------------------------------
if spm_get_defaults('stats.rft.nonstat'), rftstat = {'false','xsd:boolean'};
else                                      rftstat = {'true','xsd:boolean'}; end
idSearchSpace = getid('niiri:search_space_id',isHumanReadable);
p.entity(idSearchSpace,{...
    'prov:type',nidm_conv('nidm:SearchSpaceMaskMap'),...
    'prov:location',{uri(spm_file(files.searchspace,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.searchspace,'filename'),'xsd:string'},...
    'dct:format',niifmt,...
    'prov:label',{'Search Space Mask Map','xsd:string'}...
    nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
    nidm_conv('nidm:searchVolumeInVoxels'),{xSPM.S,'xsd:int'},...
    nidm_conv('nidm:searchVolumeInUnits'),{TabDat.ftr{8,2}(1),'xsd:float'},...
    nidm_conv('nidm:reselSizeInVoxels'),{TabDat.ftr{9,2}(end),'xsd:float'},...
    nidm_conv('nidm:searchVolumeInResels'),{TabDat.ftr{8,2}(3),'xsd:float'},...
    nidm_conv('nidm:searchVolumeReselsGeometry'),{xSPM.R,'xsd:string'},...
    nidm_conv('nidm:noiseFWHMInVoxels'),{xSPM.FWHM,'xsd:string'},...
    nidm_conv('nidm:noiseFWHMInUnits'),{TabDat.ftr{7,2}(1:3),'xsd:string'},...
    nidm_conv('nidm:randomFieldStationarity'),rftstat,...
    nidm_conv('nidm:expectedNumberOfVoxelsPerCluster'),{TabDat.ftr{3,2},'xsd:float'},...
    nidm_conv('nidm:expectedNumberOfClusters'),{TabDat.ftr{4,2},'xsd:float'},...
    nidm_conv('spm:heightCriticalThresholdFWE05'),{xSPM.uc(1),'xsd:float'},...
    nidm_conv('spm:heightCriticalThresholdFDR05'),{xSPM.uc(2),'xsd:float'},...
    nidm_conv('spm:smallestSignifClusterSizeInVoxelsFWE05'),{xSPM.uc(3),'xsd:int'},...
    nidm_conv('spm:smallestSignifClusterSizeInVoxelsFDR05'),{xSPM.uc(4),'xsd:int'},...
    'crypto:sha512',{sha512sum(spm_file(files.searchspace,'cpath')),'xsd:string'},...
    });
p.wasGeneratedBy(idSearchSpace, idInference);

%-Entity: Excursion Set
%--------------------------------------------------------------------------
if size(TabDat.dat,1) > 0
    c  = TabDat.dat{1,2};
    pc = TabDat.dat{1,1};
else
    c  = 0;
    pc = NaN;
end
idExcursionSet = getid('niiri:excursion_set_map_id',isHumanReadable);
idClusterLabelsMap = getid('niiri:cluster_label_map_id',isHumanReadable);
idMaximumIntensityProjection = getid('niiri:maximum_intensity_projection_id',isHumanReadable);
p.entity(idExcursionSet,{...
    'prov:type',nidm_conv('nidm:ExcursionSetMap'),...
    'prov:location',{uri(spm_file(files.tspm,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.tspm,'filename'),'xsd:string'},...
    'dct:format',niifmt,...
    'prov:label',{'Excursion Set Map','xsd:string'},...
    nidm_conv('nidm:numberOfClusters'),{c,'xsd:int'},...
    nidm_conv('nidm:pValue'),{pc,'xsd:float'},...
    nidm_conv('nidm:hasClusterLabelsMap'),idClusterLabelsMap,...
    nidm_conv('nidm:hasMaximumIntensityProjection'),idMaximumIntensityProjection,...
    nidm_conv('nidm:inCoordinateSpace'),id_data_coordspace,...
    'crypto:sha512',{sha512sum(spm_file(files.tspm,'cpath')),'xsd:string'},...
    });
p.entity(idClusterLabelsMap,{...
    'prov:type',nidm_conv('nidm:ClusterLabelsMap'),...
    'prov:location',{uri(spm_file(files.clust,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.clust,'filename'),'xsd:string'},...
    'dct:format',niifmt,...
    });
p.entity(idMaximumIntensityProjection,{...
    'prov:type','dctype:Image',...
    'prov:location',{uri(spm_file(files.mip,'cpath')),'xsd:anyURI'},...
    'nfo:fileName',{spm_file(files.mip,'filename'),'xsd:string'},...
    'dct:format',{'image/png','xsd:string'}...
    });
p.wasGeneratedBy(idExcursionSet, idInference); 

%-Entity: Clusters
%--------------------------------------------------------------------------
idx = find(~cellfun(@isempty,{TabDat.dat{:,5}}));
idCluster = cell(1,numel(idx));
for i=1:numel(idx)
    iClus = sprintf('%04d',i);
    idCluster{i} = getid(['niiri:significant_cluster_' iClus],isHumanReadable);
    p.entity(idCluster{i},{...
        'prov:type',nidm_conv('nidm:SignificantCluster'),...
        'prov:label',{['Significant Cluster: ' iClus],'xsd:string'},...
        nidm_conv('nidm:clusterSizeInVoxels'),{TabDat.dat{idx(i),5},'xsd:int'},...
        nidm_conv('nidm:clusterSizeInResels'),{TabDat.dat{idx(i),5}*V2R,'xsd:float'},...
        nidm_conv('nidm:pValueUncorrected'),{TabDat.dat{idx(i),6},'xsd:float'},...
        nidm_conv('nidm:pValueFWER'),{TabDat.dat{idx(i),3},'xsd:float'},...
        nidm_conv('nidm:qValueFDR'),{TabDat.dat{idx(i),4},'xsd:float'},...
        nidm_conv('nidm:clusterLabelId'),{num2str(i),'xsd:int'},...
        });
    p.wasDerivedFrom(idCluster{i}, idExcursionSet);
end

%-Entity: Peaks
%--------------------------------------------------------------------------
idx = cumsum(~cellfun(@isempty,{TabDat.dat{:,5}}));
for i=1:size(TabDat.dat,1)
    iPeak  = sprintf('%04d',i);
    idPeak = getid(['niiri:peak_' iPeak],isHumanReadable);
    idCoordinate = getid(['niiri:coordinate_' iPeak],isHumanReadable);
    
    p.entity(idPeak,{...
        'prov:type',nidm_conv('nidm:Peak'),...
        'prov:label',{['Peak: ' iPeak],'xsd:string'},...
        'prov:location',idCoordinate,...
        'prov:value',{TabDat.dat{i,9},'xsd:float'},...
        nidm_conv('nidm:equivalentZStatistic'),{xsdfloat(TabDat.dat{i,10}),'xsd:float'},...
        nidm_conv('nidm:pValueUncorrected'),{TabDat.dat{i,11},'xsd:float'},...
        nidm_conv('nidm:pValueFWER'),{TabDat.dat{i,7},'xsd:float'},...
        nidm_conv('nidm:qValueFDR'),{TabDat.dat{i,8},'xsd:float'},...
        });

    p.entity(idCoordinate,{...
        'prov:type','prov:Location',...
        'prov:type',nidm_conv('nidm:Coordinate'),...
        'prov:label',{['Coordinate: ' iPeak],'xsd:string'},...
        nidm_conv('nidm:coordinateVector'),{TabDat.dat{i,12}(1:3),'xsd:string'},...
        });
    
    p.wasDerivedFrom(idPeak, idCluster{idx(i)});
end

pp.bundle(idResults,p);

%==========================================================================
%-                  P R O V   S E R I A L I Z A T I O N
%==========================================================================
serialize(pp,fullfile(outdir,'nidm.provn'));
serialize(pp,fullfile(outdir,'nidm.ttl'));
%serialize(pp,fullfile(outdir,'nidm.json'));
%serialize(pp,fullfile(outdir,'nidm.pdf'));
%zip(fullfile(SPM.swd,[spm_file(outdir,'basename'),'.nidm.zip']),outdir)

prov = pp;


%==========================================================================
% function v = xsdfloat(v)
%==========================================================================
function v = xsdfloat(v)
% See http://books.xmlschemata.org/relaxng/ch19-77095.html
if numel(v) == 1 && isinf(v) && v > 0, v = 'INF';  end
if numel(v) == 1 && isinf(v) && v < 0, v = '-INF'; end
if numel(v) == 1 && isnan(v),          v = 'NaN';  end


%==========================================================================
% function str = esc(str)
%==========================================================================
function str = esc(str)
%-Escape
% See http://www.w3.org/TR/html4/charset.html#h-5.3.2
str = strrep(str,'&','&amp;');
str = strrep(str,'<','&lt;');
str = strrep(str,'>','&gt;');
str = strrep(str,'"','&quot;');


%==========================================================================
% function u = uri(u)
%==========================================================================
function u = uri(u)
%-File URI scheme
if ispc, s='/'; else s=''; end
u = ['file://' s strrep(spm_file(u,'cpath'),'\','/')];
e = ' ';
for i=1:length(e)
    u = strrep(u,e(i),['%' dec2hex(e(i))]);
end
u = ['file://./' spm_file(u,'filename')];


%==========================================================================
% function checksum = sha512sum(file)
%==========================================================================
function checksum = sha512sum(file)
%checksum = 'e43b6e01b0463fe7d40782137867a...'; return % TEMP (short)
md   = java.security.MessageDigest.getInstance('SHA-512');
file = spm_file(file,'cpath');
fid  = fopen(file,'rb');
if fid == -1, error('Cannot open "%s".',file); end
md.update(fread(fid,Inf,'*uint8'));
fclose(fid);
checksum = typecast(md.digest,'uint8');
checksum = lower(reshape(dec2hex(checksum)',1,[]));


%==========================================================================
% function checksum = md5sum(data)
%==========================================================================
function checksum = md5sum(data)
%checksum = 'd41d8cd98f00b204e9800998ecf8427e'; return % TEMP (short)
if ~nargin
    data = char(java.util.UUID.randomUUID);
end
md   = java.security.MessageDigest.getInstance('MD5');
if ischar(data)
    md.update(uint8(data));
else
    md.update(typecast(data,'uint8'));
end
checksum = typecast(md.digest,'uint8');
checksum = lower(reshape(dec2hex(checksum)',1,[]));


%==========================================================================
% function img2nii(img,nii,xSPM)
%==========================================================================
function img2nii(img,nii,xSPM)
if nargin == 2, xSPM = struct; end
if ~isfield(xSPM,'STAT'), xSPM.STAT = ''; end
if ~isfield(xSPM,'fcn'), xSPM.fcn = @(x) x; end
if nargin == 1, nii = spm_file(img,'ext','.nii'); end
gz = strcmp(spm_file(nii,'ext'),'gz');
if gz, nii = spm_file(nii,'ext',''); end
ni     = nifti(img);
no     = nifti;
no.dat = file_array(nii,...
                    ni.dat.dim,...
                    ni.dat.dtype,...
                    0,...
                    ni.dat.scl_slope,...
                    ni.dat.scl_inter);
no.mat  = ni.mat;
no.mat_intent = ni.mat_intent;
no.mat0 = ni.mat0;
no.mat0_intent = ni.mat0_intent;
no.descrip = ni.descrip;
switch xSPM.STAT
    case 'T'
        no.intent.name  = ['spm' xSPM.STATstr];
        no.intent.code  = 3;
        no.intent.param = xSPM.df(2);
    case 'F'
        no.intent.name  = ['spm' xSPM.STATstr];
        no.intent.code  = 4;
        no.intent.param = xSPM.df;
    case 'con'
        no.intent.name  = 'SPM contrast';
        no.intent.code  = 1001;
end

create(no);
no.dat(:,:,:) = xSPM.fcn(ni.dat(:,:,:));
if gz
    gzip(nii);
    spm_unlink(nii);
end


%==========================================================================
% function make_ROI(fname,DIM,M,xY)
%==========================================================================
function make_ROI(fname,DIM,M,xY)
gz = strcmp(spm_file(fname,'ext'),'gz');
if gz, fname = spm_file(fname,'ext',''); end
R = struct(...
    'fname',  fname,...
    'dim',    DIM,...
    'dt',     [spm_type('uint8'), spm_platform('bigend')],...
    'mat',    M,...
    'pinfo',  [1,0,0]',...
    'descrip','ROI');
Q    = zeros(DIM);
[xY, XYZmm, j] = spm_ROI(xY, struct('dim',DIM,'mat',M));
Q(j) = 1;
R    = spm_write_vol(R,Q);
if gz
    gzip(R.fname);
    spm_unlink(R.fname);
end


%==========================================================================
% function id = coordspace(p,M,DIM,units,coordsys,idx)
%==========================================================================
function id = coordspace(p,M,DIM,units,coordsys,idx)
persistent index
if nargin == 6
    index = idx;
else
    if isempty(index)
        index = 1;
    else
        index = index + 1;
    end
end
% Convert from first voxel at [1,1,1] to first voxel at [0,0,0]
v2wm = M * [eye(4,3) [1 1 1 1]'];
M    = M(1:3,1:3);
id = getid(['niiri:coordinate_space_id_' num2str(index)],isHumanReadable);
p.entity(id,{...
    'prov:type',nidm_conv('nidm:CoordinateSpace'),...
    'prov:label',{['Coordinate space ' num2str(index)],'xsd:string'},...
    nidm_conv('nidm:voxelToWorldMapping'),{v2wm,'xsd:string'},...
    nidm_conv('nidm:voxelUnits'),{units,'xsd:string'},...
    nidm_conv('nidm:voxelSize'),{sqrt(diag(M'*M))','xsd:string'},...
    nidm_conv('nidm:inWorldCoordinateSystem'),coordsys,...
    nidm_conv('nidm:numberOfDimensions'),{numel(DIM),'xsd:int'},...
    nidm_conv('nidm:dimensionsInVoxels'),{DIM,'xsd:string'}
    });

%==========================================================================
% function id = originalfile(p,file,idx,typ)
%==========================================================================
function id = originalfile(p,file,idx,typ)
id = getid([idx '_der'],isHumanReadable);
p.entity(id,{...
    'prov:type',typ,...
    'nfo:fileName',{spm_file(file,'filename'),'xsd:string'},...
    'dct:format',{'image/nifti','xsd:string'},...
    'crypto:sha512',{sha512sum(spm_file(file,'cpath')),'xsd:string'},...
    });

%==========================================================================
% function id = getid(id,humanReadable,checksum)
%==========================================================================
function id = getid(id,humanReadable,checksum)
if ~humanReadable
    if nargin == 2
        id = md5sum;
    else
        id = md5sum(checksum);
    end
end

%==========================================================================
% function i = isHumanReadable(i)
%==========================================================================
function i = isHumanReadable(i)
persistent isHR
if nargin, isHR = i; end
if isempty(isHR), error('Default not set.'); end
i = isHR;

%==========================================================================
% function 
%==========================================================================
function out = nidm_conv(in)
persistent C
if isempty(C), C = nidm_constants; end

i = find(ismember(C(:,1),in));
if ~isempty(i)
    out = C{i,2};
else
    warning('Unknown element ''%s''.',in);
    out = in;
end

%==========================================================================
% function C = nidm_constants
%==========================================================================
function C = nidm_constants
% from https://github.com/incf-nidash/nidm/blob/master/scripts/Constants.py
C = {...
'neurolex:SPM', 'neurolex:nif-0000-00343';...
'spm:partialConjunctionDegree','spm:SPM_0000015';...
'spm:smallestSignifClusterSizeInVoxelsFWE05', 'spm:SPM_0000014';...
'spm:smallestSignifClusterSizeInVoxelsFDR05', 'spm:SPM_0000013';...
'spm:smallestSignifClusterSizeInVerticessFWE05', 'spm:SPM_0000012';...
'spm:smallestSignifClusterSizeInVerticesFDR05', 'spm:SPM_0000011';...
'nidm:searchVolumeReselsGeometry', 'spm:SPM_0000010';...
'spm:TemporalDerivative', 'spm:SPM_0000006';...
'SPM_KCONJUNCTION_INFERENCE', 'spm:SPM_0000005';...
'spm:GammaDifferenceHRF', 'spm:SPM_0000004';...
'spm:DispersionDerivative', 'spm:SPM_0000003';...
'spm:DCTDriftModel', 'spm:SPM_0000002';...
'spm:driftCutoffPeriod', 'spm:SPM_0000001';...
'NIDM_CONTRAST_EXPLAINED_MEAN_SQUARE_MAP','nidm:NIDM_0000163';...
'NIDM_THRESHOLD','nidm:NIDM_0000162';...
'nidm:equivalentThreshold','nidm:NIDM_0000161';...
'nidm:pValueUncorrectedClass','nidm:NIDM_0000160';...
'nidm:noiseFWHMInVoxels', 'nidm:NIDM_0000159';...
'nidm:noiseFWHMInVertices', 'nidm:NIDM_0000158';...
'nidm:noiseFWHMInUnits', 'nidm:NIDM_0000157';...
'nidm:clusterSizeInResels','nidm:NIDM_0000156';...
'nidm:fMRIDesign','nidm:NIDM_0000155';...
'NIDM_MIXED_DESIGN','nidm:NIDM_0000154';...
'NIDM_EVENT_RELATED_DESIGN','nidm:NIDM_0000153';...
'nidm:BlockDesign','nidm:NIDM_0000152';...
'NIDM_SINE_BASIS_SET','nidm:NIDM_0000151';...
'NIDM_LINEAR_SPLINE_BASIS_SET','nidm:NIDM_0000150';...
'nidm:searchVolumeInResels','nidm:NIDM_0000149';...
'nidm:reselSizeInVoxels','nidm:NIDM_0000148';...
'spm:heightCriticalThresholdFWE05','nidm:NIDM_0000147';...
'spm:heightCriticalThresholdFDR05','nidm:NIDM_0000146';...
'NIDM_NOISE_ROUGHNESS_IN_VOXELS','nidm:NIDM_0000145';...
'nidm:ReselsPerVoxelMap','nidm:NIDM_0000144';...
'nidm:expectedNumberOfVoxelsPerCluster','nidm:NIDM_0000143';...
'nidm:expectedNumberOfVerticesPerCluster','nidm:NIDM_0000142';...
'nidm:expectedNumberOfClusters','nidm:NIDM_0000141';...
'NIDM_CLUSTER_CENTER_OF_GRAVITY','nidm:NIDM_0000140';...
'NIDM_COORDINATE_VECTOR_IN_VOXELS','nidm:NIDM_0000139';...
'nidm:hasMaximumIntensityProjection','nidm:NIDM_0000138';...
'nidm:searchVolumeInVertices','nidm:NIDM_0000137';...
'nidm:searchVolumeInUnits','nidm:NIDM_0000136';...
'NIDM_CONTRAST_VARIANCE_MAP','nidm:NIDM_0000135';...
'nidm:withEstimationMethod','nidm:NIDM_0000134';...
'nidm:voxelUnits','nidm:NIDM_0000133';...
'nidm:voxelToWorldMapping','nidm:NIDM_0000132';...
'nidm:voxelSize','nidm:NIDM_0000131';...
'nidm:voxel6connected','nidm:NIDM_0000130';...
'nidm:voxel26connected','nidm:NIDM_0000129';...
'nidm:voxel18connected','nidm:NIDM_0000128';...
'nidm:version','nidm:NIDM_0000127';...
'nidm:varianceSpatialModel','nidm:NIDM_0000126';...
'nidm:userSpecifiedThresholdType','nidm:NIDM_0000125';...
'nidm:targetIntensity','nidm:NIDM_0000124';...
'nidm:statisticType','nidm:NIDM_0000123';...
'nidm:softwareVersion','nidm:NIDM_0000122';...
'nidm:searchVolumeInVoxels','nidm:NIDM_0000121';...
'nidm:randomFieldStationarity','nidm:NIDM_0000120';...
'nidm:qValueFDR','nidm:NIDM_0000119';...
'NIDM_PIXEL8CONNECTED','nidm:NIDM_0000118';...
'NIDM_PIXEL4CONNECTED','nidm:NIDM_0000117';...
'nidm:pValueUncorrected','nidm:NIDM_0000116';...
'nidm:pValueFWER','nidm:NIDM_0000115';...
'nidm:pValue','nidm:NIDM_0000114';...
'NIDM_OBJECT_MODEL','nidm:NIDM_0000113';...
'nidm:numberOfDimensions','nidm:NIDM_0000112';...
'nidm:numberOfClusters','nidm:NIDM_0000111';...
'NIDM_GAUSSIAN_HRF','nidm:NIDM_0000110';...
'nidm:minDistanceBetweenPeaks','nidm:NIDM_0000109';...
'nidm:maxNumberOfPeaksPerCluster','nidm:NIDM_0000108';...
'nidm:maskedMedian','nidm:NIDM_0000107';...
'nidm:isUserDefined','nidm:NIDM_0000106';...
'nidm:inWorldCoordinateSystem','nidm:NIDM_0000105';...
'nidm:inCoordinateSpace','nidm:NIDM_0000104';...
'NIDM_HAS_MAP_HEADER','nidm:NIDM_0000103';...
'nidm:hasHRFBasis','nidm:NIDM_0000102';...
'nidm:hasErrorDistribution','nidm:NIDM_0000101';...
'nidm:hasErrorDependence','nidm:NIDM_0000100';...
'nidm:hasConnectivityCriterion','nidm:NIDM_0000099';...
'nidm:hasClusterLabelsMap','nidm:NIDM_0000098';...
'nidm:hasAlternativeHypothesis','nidm:NIDM_0000097';...
'nidm:grandMeanScaling','nidm:NIDM_0000096';...
'nidm:errorVarianceHomogeneous','nidm:NIDM_0000094';...
'nidm:errorDegreesOfFreedom','nidm:NIDM_0000093';...
'nidm:equivalentZStatistic','nidm:NIDM_0000092';...
'nidm:effectDegreesOfFreedom','nidm:NIDM_0000091';...
'nidm:dimensionsInVoxels','nidm:NIDM_0000090';...
'nidm:dependenceSpatialModel','nidm:NIDM_0000089';...
'nidm:hasDriftModel','nidm:NIDM_0000088';...
'NIDM_DRIFT_MODEL','nidm:NIDM_0000087';...
'nidm:coordinateVector','nidm:NIDM_0000086';...
'nidm:contrastName','nidm:NIDM_0000085';...
'nidm:clusterSizeInVoxels','nidm:NIDM_0000084';...
'nidm:clusterSizeInVertices','nidm:NIDM_0000083';...
'nidm:clusterLabelId','nidm:NIDM_0000082';...
'NIDM_WORLD_COORDINATE_SYSTEM','nidm:NIDM_0000081';...
'NIDM_VOXEL_CONNECTIVITY_CRITERION','nidm:NIDM_0000080';...
'NIDM_TWO_TAILED_TEST','nidm:NIDM_0000079';...
'NIDM_TALAIRACH_COORDINATE_SYSTEM','nidm:NIDM_0000078';...
'NIDM_SUBJECT_COORDINATE_SYSTEM','nidm:NIDM_0000077';...
'nidm:StatisticMap','nidm:NIDM_0000076';...
'NIDM_STANDARDIZED_COORDINATE_SYSTEM','nidm:NIDM_0000075';...
'nidm:SpatiallyRegularizedModel','nidm:NIDM_0000074';...
'nidm:SpatiallyLocalModel','nidm:NIDM_0000073';...
'nidm:SpatiallyGlobalModel','nidm:NIDM_0000072';...
'NIDM_SPATIAL_MODEL','nidm:NIDM_0000071';...
'nidm:SignificantCluster','nidm:NIDM_0000070';...
'NIDM_FOURIER_BASIS_SET','nidm:NIDM_0000069';...
'nidm:SearchSpaceMaskMap','nidm:NIDM_0000068';...
'NIDM_CUSTOM_BASIS_SET','nidm:NIDM_0000067';...
'nidm:ResidualMeanSquaresMap','nidm:NIDM_0000066';...
'NIDM_POISSON_DISTRIBUTION','nidm:NIDM_0000065';...
'NIDM_PIXEL_CONNECTIVITY_CRITERION','nidm:NIDM_0000064';...
'nidm:PeakDefinitionCriteria','nidm:NIDM_0000063';...
'nidm:Peak','nidm:NIDM_0000062';...
'nidm:ParameterEstimateMap','nidm:NIDM_0000061';...
'nidm:OneTailedTest','nidm:NIDM_0000060';...
'NIDM_NON_PARAMETRIC_SYMMETRIC_DISTRIBUTION','nidm:NIDM_0000059';...
'NIDM_NON_PARAMETRIC_DISTRIBUTION','nidm:NIDM_0000058';...
'nidm:objectModel','nidm:NIDM_0000057';...
'nidm:ModelParametersEstimation','nidm:NIDM_0000056';...
'NIDM_MNI305_COORDINATE_SYSTEM','nidm:NIDM_0000055';...
'nidm:MaskMap','nidm:NIDM_0000054';...
'NIDM_MAP_HEADER','nidm:NIDM_0000053';...
'NIDM_MAP','nidm:NIDM_0000052';...
'nidm:MNICoordinateSystem','nidm:NIDM_0000051';...
'NIDM_IXI549_COORDINATE_SYSTEM','nidm:NIDM_0000050';...
'nidm:Inference','nidm:NIDM_0000049';...
'nidm:IndependentError','nidm:NIDM_0000048';...
'NIDM_ICBM_MNI152_NON_LINEAR6TH_GENERATION_COORDINATE_SYSTEM','nidm:NIDM_0000047';...
'NIDM_ICBM_MNI152_NON_LINEAR2009C_SYMMETRIC_COORDINATE_SYSTEM','nidm:NIDM_0000046';...
'NIDM_ICBM_MNI152_NON_LINEAR2009C_ASYMMETRIC_COORDINATE_SYSTEM','nidm:NIDM_0000045';...
'NIDM_ICBM_MNI152_NON_LINEAR2009B_SYMMETRIC_COORDINATE_SYSTEM','nidm:NIDM_0000044';...
'NIDM_ICBM_MNI152_NON_LINEAR2009B_ASYMMETRIC_COORDINATE_SYSTEM','nidm:NIDM_0000043';...
'NIDM_ICBM_MNI152_NON_LINEAR2009A_SYMMETRIC_COORDINATE_SYSTEM','nidm:NIDM_0000042';...
'NIDM_ICBM_MNI152_NON_LINEAR2009A_ASYMMETRIC_COORDINATE_SYSTEM','nidm:NIDM_0000041';...
'NIDM_ICBM_MNI152_LINEAR_COORDINATE_SYSTEM','nidm:NIDM_0000040';...
'NIDM_ICBM452_WARP5_COORDINATE_SYSTEM','nidm:NIDM_0000039';...
'NIDM_ICBM452_AIR_COORDINATE_SYSTEM','nidm:NIDM_0000038';...
'NIDM_HEMODYNAMIC_RESPONSE_FUNCTION_DERIVATIVE','nidm:NIDM_0000037';...
'NIDM_HEMODYNAMIC_RESPONSE_FUNCTION_BASIS','nidm:NIDM_0000036';...
'NIDM_HEMODYNAMIC_RESPONSE_FUNCTION','nidm:NIDM_0000035';...
'nidm:HeightThreshold','nidm:NIDM_0000034';...
'nidm:GrandMeanMap','nidm:NIDM_0000033';...
'nidm:GaussianDistribution','nidm:NIDM_0000032';...
'NIDM_GAMMA_HRF','nidm:NIDM_0000031';...
'NIDM_GAMMA_HRB','nidm:NIDM_0000030';...
'NIDM_GAMMA_DIFFERENCE_HRF','nidm:NIDM_0000029';...
'NIDM_FINITE_IMPULSE_RESPONSE_HRB','nidm:NIDM_0000028';...
'nidm:Results','nidm:NIDM_0000027';...
'nidm:ExtentThreshold','nidm:NIDM_0000026';...
'nidm:ExcursionSetMap','nidm:NIDM_0000025';...
'NIDM_EXCHANGEABLE_ERROR','nidm:NIDM_0000024';...
'nidm:ErrorModel','nidm:NIDM_0000023';...
'NIDM_ERROR_DISTRIBUTION','nidm:NIDM_0000022';...
'nidm:regressorNames','nidm:NIDM_0000021';...
'nidm:DisplayMaskMap','nidm:NIDM_0000020';...
'nidm:DesignMatrix','nidm:NIDM_0000019';...
'nidm:Data','nidm:NIDM_0000018';...
'NIDM_CUSTOM_COORDINATE_SYSTEM','nidm:NIDM_0000017';...
'nidm:CoordinateSpace','nidm:NIDM_0000016';...
'nidm:Coordinate','nidm:NIDM_0000015';...
'NIDM_LEGENDRE_POLYNOMIAL_ORDER','nidm:NIDM_0000014';...
'nidm:ContrastStandardErrorMap','nidm:NIDM_0000013';...
'NIDM_CONNECTIVITY_CRITERION','nidm:NIDM_0000012';...
'nidm:ConjunctionInference','nidm:NIDM_0000011';...
'nidm:hasfMRIDesign','nidm:NIDM_0000010';...
'NIDM_COLIN27_COORDINATE_SYSTEM','nidm:NIDM_0000009';...
'nidm:ClusterLabelsMap','nidm:NIDM_0000008';...
'nidm:ClusterDefinitionCriteria','nidm:NIDM_0000007';...
'NIDM_CLUSTER','nidm:NIDM_0000006';...
'NIDM_BINOMIAL_DISTRIBUTION','nidm:NIDM_0000005';...
'NIDM_BINARY_MAP','nidm:NIDM_0000004';...
'NIDM_ARBITRARILY_CORRELATED_ERROR','nidm:NIDM_0000003';...
'nidm:ContrastEstimation','nidm:NIDM_0000001';...
'nidm:ContrastMap','nidm:NIDM_0000002';...
'obo:Statistic','obo:STATO_0000039';...
'obo:pValueFWER','obo:OBI_0001265';...
'OBO_Q_VALUE_FDR','obo:OBI_0001442';...
'stato:OLS','obo:STATO_0000370';...
'stato:GLS','obo:STATO_0000372';...
'stato:TStatistic','obo:STATO_0000176';...
'stato:FStatistic','obo:STATO_0000282';...
'stato:ZStatistic','obo:STATO_0000376';...
'stato:ContrastWeightMatrix','obo:STATO_0000323';...
'stato:ToeplitzCovarianceStructure','obo:STATO_0000357'};
