Skip to content
Permalink
Browse files

Merge branch 'master' into master

  • Loading branch information...
tiborauer committed Jan 31, 2019
2 parents 318ab8a + ee48221 commit c4471fc86ded29b9eb820793520848a0d45a07b3
@@ -81,10 +81,11 @@
k=1;

while (k<=size(dicomdata_subdir,1))

oldAcquisitionNumber=-1;
thispass_numvolumes=0;
DICOMHEADERS=[];

% This array is used to collect sliceing timing info so we can
% reconstruct the slice order
sliceInfoD = zeros(0, 3);
@@ -94,11 +95,11 @@
tmp = spm_dicom_headers(deblank(dicomdata_subdir(k,:)));
infoD = tmp{1};
headerFields = fieldnames(infoD);

if k == 1
chunksize_volumes=memLimit*1024/(infoD.StartOfPixelData+infoD.SizeOfPixelData)/4;
end

TR = [];
TE = [];
sliceorder = '';
@@ -137,18 +138,33 @@
end
end

% if we didn't find that private field, use standard fields.
% Philips
if isempty(TR) && isfield(infoD,'Private_2005_1030')
TR = infoD.Private_2005_1030(1);
end
if isempty(TE) && isfield(infoD,'Private_2001_1025')
TE = cellfun(@(x) str2double(x), regexp(infoD.Private_2001_1025,'[0-9.]*','match'));
end

% if we didn't find a private field, try standard fields.
if isempty(TR) && isfield(infoD, 'RepetitionTime')
TR = infoD.RepetitionTime;
else
aas_log(aap,false,'WARNING: TR not found');
end

if isempty(TE) && isfield(infoD, 'EchoTime')
TE = infoD.EchoTime;
else
aas_log(aap,false,'WARNING: TE not found');
end

% we have an obligation to warn the user
if isempty(TR)
aas_log(aap,false,'WARNING: TR not found');
TR = NaN;
end
if isempty(TE)
aas_log(aap,false,'WARNING: TE not found, setting to NaN');
TE = NaN;
end

% Siemens
if isempty(sliceorder) && isfield(infoD, 'CSAImageHeaderInfo') && cell_index({infoD.CSAImageHeaderInfo.name},'MosaicRefAcqTimes')
slicetimes = aas_get_numaris4_numval(infoD.CSAImageHeaderInfo,'MosaicRefAcqTimes')';
@@ -178,7 +194,7 @@
if ~any(strcmpi(headerFields,'NumberOfPhaseEncodingSteps')), infoD.NumberOfPhaseEncodingSteps = infoD.AcquisitionMatrix(1); end
collectSOinfo = isempty(sliceorder) && isfield(infoD, 'TemporalPositionIdentifier');

% Philips (based on Hester Breman's code for BrainVoyager)
% Philips (based on Hester Breman's code for BrainVoyager)
if isempty(echospacing) && all(isfield(infoD,{'EPIFactor' 'WaterFatShift' 'MagneticFieldStrength'}))
epifactor = infoD.EPIFactor;
water_fat_shift_pixel = infoD.WaterFatShift;
@@ -200,13 +216,32 @@
infoD.slicetimes = slicetimes/1000;
infoD.echospacing = echospacing;

% Single slice per DICOM:
% TemporalPositionIdentifier is basically the volume number
% InstanceNumber is the temporal position in that acqusition
% SliceLocation is spatial
sliceInfoD(end+1, 2:3) = [infoD.InstanceNumber infoD.SliceLocation];
if collectSOinfo
sliceInfoD(end, 1) = infoD.TemporalPositionIdentifier;
% sanity check: some aa modules assume fieldnames
% RepetitionTime and EchoTime exist in the DICOM header
% which might not if they were defined by private
% fields -- create them now if need be (note also
% these are in milliseconds not seconds)

if (~isfield(infoD,'EchoTime'))
infoD.EchoTime = TE;
end

if (~isfield(infoD,'RepetitionTime'))
infoD.RepetitionTime = TR;
end

if isfield(infoD,'SliceLocation')

% Single slice per DICOM:
% TemporalPositionIdentifier is basically the volume number
% InstanceNumber is the temporal position in that acqusition
% SliceLocation is spatial

sliceInfoD(end+1, 2:3) = [infoD.InstanceNumber infoD.SliceLocation];
if collectSOinfo
sliceInfoD(end, 1) = infoD.TemporalPositionIdentifier;
end

end

DICOMHEADERS=[DICOMHEADERS {infoD}];
@@ -244,9 +279,9 @@
end

aas_log(aap,false,sprintf('Sliceorder %s have been detected', sliceorder));

% Update the DICOMHEADERS
[junk, sliceorder] = sort(sliceInfo(:,2)); % for Philips
[junk, sliceorder] = sort(sliceInfo(:,2)); % for Philips
DICOMHEADERS = arrayfun(@(x) {setfield(x{1}, 'sliceorder', sliceorder')}, DICOMHEADERS);
DICOMHEADERS = arrayfun(@(x) {setfield(x{1}, 'slicetimes', x{1}.volumeTR/numSlices*(x{1}.sliceorder-1))}, DICOMHEADERS);
end
@@ -261,6 +296,7 @@
end
end
end

% [AVG] to cope with modern cutting edge scanners, and other probs
% (e.g. 7T Siemens scanners, which seem to mess up the ICE dimensions...)
if isfield(aap.directory_conventions, 'dicom_converter') && ~isempty(aap.directory_conventions.dicom_converter)
@@ -277,36 +313,43 @@
dicom_converter = 'spm_dicom_convert';
opts = {'all','flat','nii'};
end

conv = feval(dicom_converter,DICOMHEADERS_selected,opts{:});
if ~isempty(custompath), rmpath(custompath); end
out=[out(:);conv.files(:)];

% one DICOM per slice --> one header per volume
firstsliceInd = sliceInfoD(:,3) == sliceInfoD(find(sliceInfoD(:,2)==min(sliceInfoD(:,2)),1,'first'),3); % select the ones for the first slices per volume
DICOMHEADERS = DICOMHEADERS(firstsliceInd);
[junk, sortind] = sort(sliceInfoD(firstsliceInd,2));

if ~isempty(sliceInfoD)
firstsliceInd = sliceInfoD(:,3) == sliceInfoD(find(sliceInfoD(:,2)==min(sliceInfoD(:,2)),1,'first'),3); % select the ones for the first slices per volume
DICOMHEADERS = DICOMHEADERS(firstsliceInd);
[junk, sortind] = sort(sliceInfoD(firstsliceInd,2));
else
sortind = 1;
end

DICOMHEADERS = DICOMHEADERS(sortind);

dicomheader{subdirind}=[dicomheader{subdirind} DICOMHEADERS];
end
out_allechoes{subdirind}=unique(out);

% if ~isempty(strfind(inputstream, 'dicom_structural'))
% % [AVG] This is to cope with a number of strucutral images, so we
% % may have the DICOM header of each of them...
% InstanceNumbers = [];
% DCMnumbers = [];
% % Loop throught the dicoms to see if the SeriesNumber changes
% for l=1:length(DICOMHEADERS);
% InstanceNumber = 100*DICOMHEADERS{l}.SeriesNumber + DICOMHEADERS{l}.EchoNumbers;
% if all(InstanceNumbers ~= InstanceNumber)
% InstanceNumbers(end+1) = InstanceNumber;
% DCMnumbers = [DCMnumbers l];
% end
% end
% [junk, so] = sort(InstanceNumbers);
% dicomheader={DICOMHEADERS{DCMnumbers(so)}};
% end
% if ~isempty(strfind(inputstream, 'dicom_structural'))
% % [AVG] This is to cope with a number of strucutral images, so we
% % may have the DICOM header of each of them...
% InstanceNumbers = [];
% DCMnumbers = [];
% % Loop throught the dicoms to see if the SeriesNumber changes
% for l=1:length(DICOMHEADERS);
% InstanceNumber = 100*DICOMHEADERS{l}.SeriesNumber + DICOMHEADERS{l}.EchoNumbers;
% if all(InstanceNumbers ~= InstanceNumber)
% InstanceNumbers(end+1) = InstanceNumber;
% DCMnumbers = [DCMnumbers l];
% end
% end
% [junk, so] = sort(InstanceNumbers);
% dicomheader={DICOMHEADERS{DCMnumbers(so)}};
% end
end

if numel(dicomheader) == 1
@@ -22,9 +22,11 @@
end

try val = val{index};
catch
aas_log(aap,false,sprintf('WARNING: Setting <%s> has fewer then %d elements.\nWARNING: First value will be applied!',settingstring,index));
catch E
aas_log(aap,false,sprintf('WARNING (%s): %s requested %s(%d), but', mfilename, E.stack(min(2,numel(E.stack))).name, settingstring, index));
aas_log(aap,false,sprintf('WARNING (%s): only %d value(s) are defined in the current settings.', mfilename, numel(val)));
aas_log(aap,false,sprintf('WARNING (%s): The first value (%0.9g) will be returned instead.', mfilename, val{1}));
val = val{1};
end
end
end
end
@@ -15,7 +15,7 @@
resp=sprintf('Performing dicom to nifti conversion of EPIs');

case 'summary'
resp=sprintf('Made directory %s and did dicom to nifti conversion of EPIs\n');
resp=sprintf('Made directory and did dicom to nifti conversion of EPIs\n');

case 'report'

@@ -419,22 +419,35 @@
aas_makedir(aap,aas_getpath_bydomain(aap,domain,indices));

[aap fns DICOMHEADERS]=aas_convertseries_fromstream(aap,domain,indices,'dicom_epi');

if ~isempty(aas_getsetting(aap,'ignoreafter',sess))
fns = fns(1:min(numel(fns),aas_getsetting(aap,'ignoreafter',sess)));
end
if ~isempty(aas_getsetting(aap,'ignoreafter'))
fns = fns(1:min(numel(fns),aas_getsetting(aap,'ignoreafter',sess)));
end

finalepis=fns;

if aap.tasklist.currenttask.settings.outputstats || aap.options.NIFTI4D

% Find temporal SNR
for fileind=1:numel(fns)
V(fileind)=spm_vol(fns{fileind});
end

if (numel(fns) > 1)
for fileind=1:numel(fns)
V(fileind)=spm_vol(fns{fileind});
end
else
temp=spm_vol(fns{1});
for fileind=1:numel(temp)
V(fileind)=temp(fileind);
end
end

V0 = V; % save V for 4D conversion [TA]

end;

if aap.tasklist.currenttask.settings.outputstats
[Y XYZ]=spm_read_vols(V);

[Y XYZ]=spm_read_vols(V);
for ind=1:numel(V)
if (ind==1)
Ytot=Y(:,:,:,1);
@@ -512,13 +525,17 @@
finalepis = finalepis(d+1:end);

% 4D conversion [TA]
if isfield(aap.options, 'NIFTI4D') && aap.options.NIFTI4D

if numel(finalepis) > 1 && isfield(aap.options, 'NIFTI4D') && aap.options.NIFTI4D
finalepis = finalepis{1};
ind = find(finalepis=='-');
if numel(ind) > 1, ind = ind(2); end
finalepis = [finalepis(1:ind-1) '.nii'];
spm_file_merge(char({V0(ndummies+1:end).fname}),finalepis,0,DICOMHEADERS{1}.volumeTR);
end
else
finalepis = finalepis{1};
end

% And describe outputs
aap = aas_desc_outputs(aap,domain,indices,'dummyscans',dummylist);
dcmhdrfn = fullfile(domainpath,'dicom_headers.mat');
@@ -15,6 +15,10 @@
</echoweightmethod>

<outputstats>1</outputstats> <!-- output noise estimates and so on? -->

<!-- ignoreafter can be a vector (one value for each session) -->
<!-- if scalar, the value will be applied to all sessions -->
<!-- (and this may generate a warning during execution) -->

<ignoreafter desc='Only take this many volume (including dummies)'>9999999</ignoreafter>

0 comments on commit c4471fc

Please sign in to comment.
You can’t perform that action at this time.