segment_claustrum_rh.m

function segment_claustrum_rh(t1wFile, tt5File, fslSegFile, outFile)
% extract_claustrum_rh segments the dorsal clausgrum of the right
% hemisphere from a high-resolution T1-weighted file. Segmentation is based
% on FSL's FIRST segmentation of the putamen, and on a mask of the CSF, as
% identified using MRTrix's 5ttgen.
% The function was developed for HCP's Young Adult dataset, with high
% resolution data (7^3 mm).
%
% INPUT
% t1wFile: T1-weighted file (e.g., SUBJECT/T1w/T1w_acpc_dc_restore.nii.gz)
% tt5File: A 5tt (five tissue type) file generated by MRTrix's 5ttgen, with
%          the T1-weighted file as input
% fslSegFile: output of FSL's FIRST (e.g., first_all_none_firstseg.nii.gz)
% outFile: Name of the output file to be generated
%
% OUTPUT
% outFile is written to disk, containing the segmented dorsal claustrum
% Additional intermediate files are also written, and can be used for
% debugging.
%
% Requirements: Vistasoft's function readFileNifti.m

%% Check that all input files exist
if ~(exist(t1wFile,'file') && (exist(tt5File,'file') || exist(fslSegFile,'file')))
    error('Subject is missing at least one of the necessary input files');
end

outDir = fileparts(outFile);
if ~exist(outDir,'dir') % Make output direcoty
    mkdir(outDir)
end


%% (1) Load input files
t1w = readFileNifti(t1wFile);
full5tt = readFileNifti(tt5File);

% Extract putamen
tmp = readFileNifti(fslSegFile);
putamen = zeros(size(tmp.data));
putamen(tmp.data==51) = 1;
tmp.data = putamen;
putamen = tmp;
clear tmp

%% (2) Find the Putamen edge and expand to get candidate ROI
NumVoxExp = round(20 / t1w.pixdim(1));
putamenExpanded = putamen.data;
for v = 1:NumVoxExp
    putamenExpanded(find(circshift(putamenExpanded ,1,1))) = 1;
end

putamenEdge = putamenExpanded;
putamenEdge(find(circshift(putamenExpanded,1,1))) = 0;
dtiWriteNiftiWrapper(putamenEdge, putamen.qto_xyz, fullfile(outDir,'R_putamenEdge.nii.gz'));

% Expand the Putamen edge to define candidate voxels
NumVoxExp = round(5 / t1w.pixdim(1));
candidate = zeros(size(putamen.data));
for v = 2:NumVoxExp
    candidate(find(circshift(putamenEdge,-v,1))) = 1;
end
dtiWriteNiftiWrapper(candidate, putamen.qto_xyz, fullfile(outDir,'R_candidate1.nii.gz'));


%% (3) Remove voxels that are too close to CSF
% We cannot remove the insula itself, since the insular segmentation is
% often corrupted and includes voxels of the claustrum or the extreme
% capsules. Instead, we use the CSF segmentation as follows.

% Claustrum indices
tmp = candidate;
[x,y,z] = ind2sub(size(tmp),find(tmp));
% CSF indices
csfMask = full5tt.data(:,:,:,4);
csfMask = single(csfMask>0.1);
% Keep only CSF indices lateral to the claustrum
NumVoxExp = round(5 / t1w.pixdim(1));
for v=1:NumVoxExp
    tmp(find(circshift(tmp,v,1))) = 1;
end
csfMask(find(tmp)) = 0;

[xcsf,ycsf,zcsf] = ind2sub(size(csfMask),find(csfMask));

D = pdist2([xcsf,ycsf,zcsf],[x,y,z],'euclidean','smallest',1);
highDistInd = find(D>5); % 5 mm away from CSF
goodInd = sub2ind(size(tmp),x(highDistInd ),y(highDistInd ),z(highDistInd ));

tmp = zeros(size(candidate));
tmp(goodInd)=1;
dtiWriteNiftiWrapper(tmp, putamen.qto_xyz, [outDir,'/R_candidate_farFromCsf.nii.gz']);
candidate = tmp;

%% (4) k-means clustering
t1wCandidate = t1w.data(find(candidate));
[kidx,kclust] = kmeans(t1wCandidate,2,'MaxIter',1000);
if kclust(1)>kclust(2) % make sure the 1st cluster has lower values
    kidx = 1./kidx*2;
    kclust = circshift(kclust,1,1);
end
candidateIndices = find(candidate);
kmeanLow = zeros(size(candidate));
kmeanLow(candidateIndices(kidx==1)) = 1;
dtiWriteNiftiWrapper(kmeanLow, putamen.qto_xyz, [outDir,'/R_clstrm1_kmeansLow.nii.gz']);
kmeanHigh = zeros(size(candidate));
kmeanHigh(candidateIndices(kidx==2)) = 2;
dtiWriteNiftiWrapper(kmeanHigh, putamen.qto_xyz, [outDir,'/R_kmeansHigh.nii.gz']);

%% (5) Smooth the Claustrum voxels
filt = 'box';  % 'gaussian', 'box'
sz = [1,3,3];
smoothedClaustrum = smooth3(kmeanLow,filt,sz);
smoothedClaustrum(smoothedClaustrum<0.25) = 0;
smoothedClaustrum = single(logical(smoothedClaustrum + kmeanLow));
dtiWriteNiftiWrapper(smoothedClaustrum, putamen.qto_xyz, [outDir,'/R_clstrm2_kmeansLow_smooth25.nii.gz']);

% Keep only voxels within the original candidate ROI
smoothedClaustrum(find(~candidate)) = 0;
dtiWriteNiftiWrapper(smoothedClaustrum, putamen.qto_xyz, [outDir,'/R_clstrm3_kmeansLow_smooth25_inCandidates.nii.gz']);


%% (6) keep only biggest clusters
CC = bwconncomp(smoothedClaustrum);
ccLen = cellfun(@length,CC.PixelIdxList);
idx = find(ccLen>20);

claustrum= zeros(size(candidate));
for ii=1:length(idx)
    claustrum(CC.PixelIdxList{idx(ii)}) = 1;
end

%% (7) Save
% Write final result to the subject's directory
dtiWriteNiftiWrapper(claustrum, putamen.qto_xyz, outFile);

end