On Wed, 7 Sep 2005, Toshihiko Aso MD wrote:
> Is there an easier way to import DICOM files into the original
> folders they are stored, instead of PWD? It looks intricate to
> modify spm_dicom_convert. I am using SPM5b.
With some Matlab programming expertise, it is not too intricate to do some
modifications to spm_dicom_convert. You would only have to find the
places where the output filenames are constructed and modify the paths and
filenames according to your convention. Personally, I would not mix up
DICOM and analyze data for at least 2 reasons:
- It is good to keep your DICOM repository "clean" - this is your raw data
which is the foundation of any analysis.
- There are so many DICOM files (especially if you are using non-mosaic
sequences) per directory, that file selection will be rather slow in a
folder with DICOM and nifti files.
I have attached a modified version of spm_dicom_convert and
spm_config_dicom.m that may give you some idea on how to implement your
own path name strategies. I have moved the code to determine the path and
filename into a separate subfunction "getfilelocation" - have a look at
this part.
Hope this helps,
--
Volkmar Glauche
-
Department of Neurology [log in to unmask]
Universitaetsklinikum Freiburg Phone 49(0)761-270-5331
Breisacher Str. 64 Fax 49(0)761-270-5416
79106 Freiburg
function spm_dicom_convert(hdr,opts,root_dir)
% Convert DICOM images into something that SPM can use
% FORMAT spm_dicom_convert(hdr,opts)
% hdr - a cell array of DICOM headers from spm_dicom_headers
% opts - options
% 'all' - all DICOM files (default)
% 'mosaic' - the mosaic images
% 'standard' - standard DICOM files
% 'raw' - convert raw FIDs (not implemented)
%
% Converted files are written to the current directory
%_______________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
% @(#)Hacked by VG to sort volumes by ICE_dims and read diffusion information
% diffusion information: flip y direction of gradient vector
% John Ashburner && Jesper Andersson
% $Id: spm_dicom_convert.m,v 1.4 2005/08/18 13:07:13 glauche Exp $
if nargin<2, opts = 'all'; end;
if nargin<3, root_dir='projects_cur';end;
[images,guff] = select_tomographic_images(hdr);
[mosaic,standard] = select_mosaic_images(images);
if (strcmp(opts,'all') || strcmp(opts,'mosaic')) && ~isempty(mosaic),
convert_mosaic(mosaic,root_dir);
end;
if (strcmp(opts,'all') || strcmp(opts,'standard')) && ~isempty(standard),
convert_standard(standard,root_dir);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function convert_mosaic(hdr,root_dir)
spm_progress_bar('Init',length(hdr),'Writing Mosaic', 'Files written');
for i=1:length(hdr),
% Output filename
%-------------------------------------------------------------------
if ~isfield(hdr{i},'ProtocolName')
if isfield(hdr{1},'SequenceName')
hdr{1}.ProtocolName = hdr{1}.SequenceName;
else
hdr{1}.ProtocolName='unknown';
end;
end;
if ~isfield(hdr{i},'EchoNumbers')
hdr{i}.EchoNumbers = 0;
end;
fname = getfilelocation(hdr{i},root_dir);
% Image dimensions and data
%-------------------------------------------------------------------
nc = hdr{i}.Columns;
nr = hdr{i}.Rows;
dim = [0 0 0];
dim(3) = read_NumberOfImagesInMosaic(hdr{i});
np = [nc nr]/ceil(sqrt(dim(3)));
dim(1:2) = np;
if ~all(np==floor(np)),
warning('%s: dimension problem [Num Images=%d, Num Cols=%d, Num Rows=%d].',...
hdr{i}.Filename,dim(3), nc,nr);
continue;
end;
% Apparently, this is not the right way of doing it.
%np = read_AcquisitionMatrixText(hdr{i});
%if rem(nc, np(1)) || rem(nr, np(2)),
% warning('%s: %dx%d wont fit into %dx%d.',hdr{i}.Filename,...
% np(1), np(2), nc,nr);
% return;
%end;
%dim = [np read_NumberOfImagesInMosaic(hdr{i})];
mosaic = read_image_data(hdr{i});
volume = zeros(dim);
for j=1:dim(3),
img = mosaic((1:np(1))+np(1)*rem(j-1,nc/np(1)), (np(2):-1:1)+np(2)*floor((j-1)/(nc/np(1))));
if ~any(img(:)),
volume = volume(:,:,1:(j-1));
break;
end;
volume(:,:,j) = img;
end;
dim = size(volume);
dt = [spm_type('int16') spm_platform('bigend')];
% Orientation information
%-------------------------------------------------------------------
% Axial Analyze voxel co-ordinate system:
% x increases right to left
% y increases posterior to anterior
% z increases inferior to superior
% DICOM patient co-ordinate system:
% x increases right to left
% y increases anterior to posterior
% z increases inferior to superior
% T&T co-ordinate system:
% x increases left to right
% y increases posterior to anterior
% z increases inferior to superior
analyze_to_dicom = [diag([1 -1 1]) [0 (dim(2)+1) 0]'; 0 0 0 1]*[eye(4,3) [-1 -1 -1 1]'];
vox = [hdr{i}.PixelSpacing hdr{i}.SpacingBetweenSlices];
pos = hdr{i}.ImagePositionPatient';
orient = reshape(hdr{i}.ImageOrientationPatient,[3 2]);
orient(:,3) = null(orient');
if det(orient)<0, orient(:,3) = -orient(:,3); end;
% The image position vector is not correct. In dicom this vector points to
% the upper left corner of the image. Perhaps it is unlucky that this is
% calculated in the syngo software from the vector pointing to the center of
% the slice (keep in mind: upper left slice) with the enlarged FoV.
dicom_to_patient = [orient*diag(vox) pos ; 0 0 0 1];
truepos = dicom_to_patient *[(size(mosaic)-dim(1:2))/2 0 1]';
dicom_to_patient = [orient*diag(vox) truepos(1:3) ; 0 0 0 1];
patient_to_tal = diag([-1 -1 1 1]);
mat = patient_to_tal*dicom_to_patient*analyze_to_dicom;
% Maybe flip the image depending on SliceNormalVector from 0029,1010
%-------------------------------------------------------------------
SliceNormalVector = read_SliceNormalVector(hdr{i});
if det([reshape(hdr{i}.ImageOrientationPatient,[3 2]) SliceNormalVector(:)])<0;
volume = volume(:,:,end:-1:1);
mat = mat*[eye(3) [0 0 -(dim(3)-1)]'; 0 0 0 1];
end;
% Possibly useful information
%-------------------------------------------------------------------
tim = datevec(hdr{i}.AcquisitionTime/(24*60*60));
descrip = sprintf('%gT %s %s TR=%gms/TE=%gms/FA=%gdeg %s %d:%d:%.5g Mosaic',...
hdr{i}.MagneticFieldStrength, hdr{i}.MRAcquisitionType,...
deblank(hdr{i}.ScanningSequence),...
hdr{i}.RepetitionTime,hdr{i}.EchoTime,hdr{i}.FlipAngle,...
datestr(hdr{i}.AcquisitionDate),tim(4),tim(5),tim(6));
% descrip = [deblank(descrip) ' ' hdr{i}.PatientsName];
if ~true, % LEFT-HANDED STORAGE
mat = mat*[-1 0 0 (dim(1)+1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
volume = flipdim(volume,1);
end;
%if isfield(hdr{i},'RescaleSlope') && hdr{i}.RescaleSlope ~= 1,
% volume = volume*hdr{i}.RescaleSlope;
%end;
%if isfield(hdr{i},'RescaleIntercept') && hdr{i}.RescaleIntercept ~= 0,
% volume = volume + hdr{i}.RescaleIntercept;
%end;
%V = struct('fname',fname, 'dim',dim, 'dt',dt, 'mat',mat, 'descrip',descrip);
%spm_write_vol(V,volume);
% Note that data are no longer scaled by the maximum amount.
% This may lead to rounding errors in smoothed data, but it
% will get around other problems.
RescaleSlope = 1;
RescaleIntercept = 0;
if isfield(hdr{i},'RescaleSlope') && hdr{i}.RescaleSlope ~= 1,
RescaleSlope = hdr{i}.RescaleSlope;
end;
if isfield(hdr{i},'RescaleIntercept') && hdr{i}.RescaleIntercept ~= 0,
RescaleIntercept = hdr{i}.RescaleIntercept;
end;
N = nifti;
N.dat = file_array(fname,dim,dt,0,RescaleSlope,RescaleIntercept);
N.mat = mat;
N.mat0 = mat;
N.mat_intent = 'Scanner';
N.mat0_intent = 'Scanner';
N.descrip = descrip;
create(N);
% Write the data unscaled
dat = file_array(fname,dim,dt,0);
for i=1:dim(3),
dat(:,:,i) = volume(:,:,i);
end;
spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
return;
%_______________________________________________________________________
%_______________________________________________________________________
function convert_standard(hdr,root_dir)
hdr = sort_into_volumes(hdr);
for i=1:length(hdr),
write_volume(hdr{i},root_dir);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function vol = sort_into_volumes(hdr)
%
% First of all, sort into volumes based on relevant
% fields in the header.
%
vol{1}{1} = hdr{1};
for i=2:length(hdr),
orient = reshape(hdr{i}.ImageOrientationPatient,[3 2]);
xy1 = hdr{i}.ImagePositionPatient*orient;
match = 0;
if isfield(hdr{i},'CSAImageHeaderInfo')
ice1 = sscanf( ...
strrep(get_numaris4_val(hdr{i}.CSAImageHeaderInfo,'ICE_Dims'), ...
'X', '-1'), '%i_%i_%i_%i_%i_%i_%i_%i_%i')';
dimsel = logical([1 1 1 1 1 1 0 0 1]);
else
ice1 = [];
end;
for j=1:length(vol),
orient = reshape(vol{j}{1}.ImageOrientationPatient,[3 2]);
xy2 = vol{j}{1}.ImagePositionPatient*orient;
dist2 = sum((xy1-xy2).^2);
if strcmp(hdr{i}.Modality,'CT') && ...
strcmp(vol{j}{1}.Modality,'CT') % Our CT seems to
% have shears in
% slice positions
dist2 = 0;
end;
if ~isempty(ice1) && isfield(vol{j}{1},'CSAImageHeaderInfo')
% Replace 'X' in ICE_Dims by '-1'
ice2 = sscanf( ...
strrep(get_numaris4_val(vol{j}{1}.CSAImageHeaderInfo,'ICE_Dims'), ...
'X', '-1'), '%i_%i_%i_%i_%i_%i_%i_%i_%i')';
if ~isempty(ice2)
identical_ice_dims=all(ice1(dimsel)==ice2(dimsel));
else
identical_ice_dims = 0; % have ice1 but not ice2, ->
% something must be different
end,
else
identical_ice_dims = 1; % No way of knowing if the is no CSAImageHeaderInfo
end;
try
match = hdr{i}.SeriesNumber == vol{j}{1}.SeriesNumber &&...
hdr{i}.Rows == vol{j}{1}.Rows &&...
hdr{i}.Columns == vol{j}{1}.Columns &&...
all(hdr{i}.ImageOrientationPatient == vol{j}{1}.ImageOrientationPatient) &&...
all(hdr{i}.PixelSpacing == vol{j}{1}.PixelSpacing) && dist2<0.00001 &&...
identical_ice_dims;
if (hdr{i}.AcquisitionNumber ~= hdr{i}.InstanceNumber)|| ...
(vol{j}{1}.AcquisitionNumber ~= vol{j}{1}.InstanceNumber)
match = match && (hdr{i}.AcquisitionNumber == vol{j}{1}.AcquisitionNumber);
end;
if isfield(hdr{i},'SequenceName') && isfield(vol{j}{1}, ...
'SequenceName')
match = match && strcmp(hdr{i}.SequenceName, ...
vol{j}{1}.SequenceName);
end;
if isfield(hdr{i},'SeriesInstanceUID') && isfield(vol{j}{1}, ...
'SeriesInstanceUID')
match = match && strcmp(hdr{i}.SeriesInstanceUID, ...
vol{j}{1}.SeriesInstanceUID);
end;
if isfield(hdr{i},'EchoNumbers') && isfield(vol{j}{1}, ...
'EchoNumbers')
match = match && hdr{i}.EchoNumbers == ...
vol{j}{1}.EchoNumbers;
end;
catch
match = 0;
end
if match
vol{j}{end+1} = hdr{i};
match = 1;
break;
end;
end;
if ~match,
vol{end+1}{1} = hdr{i};
end;
end;
%
% Secondly, sort volumes into ascending/descending
% slices depending on .ImageOrientationPatient field.
%
vol2 = {};
for j=1:length(vol),
orient = reshape(vol{j}{1}.ImageOrientationPatient,[3 2]);
proj = null(orient');
if det([orient proj])<0, proj = -proj; end;
z = zeros(length(vol{j}),1);
for i=1:length(vol{j}),
z(i) = vol{j}{i}.ImagePositionPatient*proj;
end;
[z,index] = sort(z);
vol{j} = vol{j}(index);
if length(vol{j})>1,
% dist = diff(z);
if any(diff(z)==0)
tmp = sort_into_vols_again(vol{j});
vol{j} = tmp{1};
vol2 = {vol2{:} tmp{2:end}};
end;
end;
end;
vol = {vol{:} vol2{:}};
for j=1:length(vol),
if length(vol{j})>1,
orient = reshape(vol{j}{1}.ImageOrientationPatient,[3 2]);
proj = null(orient');
if det([orient proj])<0, proj = -proj; end;
z = zeros(length(vol{j}),1);
for i=1:length(vol{j}),
z(i) = vol{j}{i}.ImagePositionPatient*proj;
end;
[z,index] = sort(z);
dist = diff(z);
if sum((dist-mean(dist)).^2)/length(dist)>0.00001,
fprintf('***************************************************\n');
fprintf('* VARIABLE SLICE SPACING *\n');
fprintf('* This may be due to missing DICOM files. *\n');
if checkfields(vol{j}{1},'PatientID','SeriesNumber','AcquisitionNumber','InstanceNumber'),
fprintf('* %s / %d / %d / %d \n',...
deblank(vol{j}{1}.PatientID), vol{j}{1}.SeriesNumber, ...
vol{j}{1}.AcquisitionNumber, vol{j}{1}.InstanceNumber);
fprintf('* *\n');
end;
fprintf('* %20.3g *\n', dist);
fprintf('***************************************************\n');
end;
end;
end;
%dcm = vol;
%save('dicom_headers.mat','dcm');
return;
%_______________________________________________________________________
%_______________________________________________________________________
function vol2 = sort_into_vols_again(volj)
if ~isfield(volj{1},'InstanceNumber'),
fprintf('***************************************************\n');
fprintf('* The slices may be all mixed up and the data *\n');
fprintf('* not really usable. Talk to your physicists *\n');
fprintf('* about this. *\n');
fprintf('***************************************************\n');
vol2 = {volj};
return;
end;
fprintf('***************************************************\n');
fprintf('* The AcquisitionNumber counter does not appear *\n');
fprintf('* to be changing from one volume to another. *\n');
fprintf('* Another possible explanation is that the same *\n');
fprintf('* DICOM slices are used multiple times. *\n');
%fprintf('* Talk to your MR sequence developers or scanner *\n');
%fprintf('* supplier to have this fixed. *\n');
fprintf('* The conversion is having to guess how slices *\n');
fprintf('* should be arranged into volumes. *\n');
if checkfields(volj{1},'PatientID','SeriesNumber','AcquisitionNumber'),
fprintf('* %s / %d / %d\n',...
deblank(volj{1}.PatientID), volj{1}.SeriesNumber, ...
volj{1}.AcquisitionNumber);
end;
fprintf('***************************************************\n');
z = zeros(length(volj),1);
t = zeros(length(volj),1);
d = zeros(length(volj),1);
orient = reshape(volj{1}.ImageOrientationPatient,[3 2]);
proj = null(orient');
if det([orient proj])<0, proj = -proj; end;
for i=1:length(volj),
z(i) = volj{i}.ImagePositionPatient*proj;
t(i) = volj{i}.InstanceNumber;
end;
% msg = 0;
[t,index] = sort(t);
volj = volj(index);
z = z(index);
msk = find(diff(t)==0);
if any(msk),
% fprintf('***************************************************\n');
% fprintf('* These files have the same InstanceNumber: *\n');
% for i=1:length(msk),
% [tmp,nam1,ext1] = fileparts(volj{msk(i)}.Filename);
% [tmp,nam2,ext2] = fileparts(volj{msk(i)+1}.Filename);
% fprintf('* %s%s = %s%s (%d)\n', nam1,ext1,nam2,ext2, volj{msk(i)}.InstanceNumber);
% end;
% fprintf('***************************************************\n');
index = [true ; diff(t)~=0];
t = t(index);
z = z(index);
d = d(index);
volj = volj(index);
end;
%if any(diff(sort(t))~=1), msg = 1; end;
[z,index] = sort(z);
volj = volj(index);
t = t(index);
vol2 = {};
while ~all(d),
i = find(~d);
i = i(1);
i = find(z==z(i));
[t(i),si] = sort(t(i));
volj(i) = volj(i(si));
for i1=1:length(i),
if length(vol2)<i1, vol2{i1} = {}; end;
vol2{i1} = {vol2{i1}{:} volj{i(i1)}};
end;
d(i) = 1;
end;
msg = 0;
if any(diff(sort(t))~=1), msg = 1; end;
if ~msg,
len = length(vol2{1});
for i=2:length(vol2),
if length(vol2{i}) ~= len,
msg = 1;
break;
end;
end;
end;
if msg,
fprintf('***************************************************\n');
fprintf('* There are missing DICOM files, so the the *\n');
fprintf('* resulting volumes may be messed up. *\n');
if checkfields(volj{1},'PatientID','SeriesNumber','AcquisitionNumber'),
fprintf('* %s / %d / %d\n',...
deblank(volj{1}.PatientID), volj{1}.SeriesNumber, ...
volj{1}.AcquisitionNumber);
end;
fprintf('***************************************************\n');
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function write_volume(hdr,root_dir)
% Output filename
%-------------------------------------------------------------------
if ~isfield(hdr{1},'ProtocolName')
if isfield(hdr{1},'SequenceName')
hdr{1}.ProtocolName = hdr{1}.SequenceName;
else
hdr{1}.ProtocolName='unknown';
end;
end;
if ~isfield(hdr{1},'EchoNumbers')
hdr{1}.EchoNumbers = 0;
end;
fname = getfilelocation(hdr{1}, root_dir);
% Image dimensions
%-------------------------------------------------------------------
nc = hdr{1}.Columns;
nr = hdr{1}.Rows;
dim = [nc nr length(hdr)];
dt = [spm_type('int16') spm_platform('bigend')];
% Orientation information
%-------------------------------------------------------------------
% Axial Analyze voxel co-ordinate system:
% x increases right to left
% y increases posterior to anterior
% z increases inferior to superior
% DICOM patient co-ordinate system:
% x increases right to left
% y increases anterior to posterior
% z increases inferior to superior
% T&T co-ordinate system:
% x increases left to right
% y increases posterior to anterior
% z increases inferior to superior
analyze_to_dicom = [diag([1 -1 1]) [0 (dim(2)+1) 0]'; 0 0 0 1]*[eye(4,3) [-1 -1 -1 1]'];
orient = reshape(hdr{1}.ImageOrientationPatient,[3 2]);
orient(:,3) = null(orient');
if det(orient)<0, orient(:,3) = -orient(:,3); end;
if length(hdr)>1,
z = zeros(length(hdr),1);
for i=1:length(hdr),
z(i) = hdr{i}.ImagePositionPatient*orient(:,3);
end;
z = mean(diff(z));
else
if checkfields(hdr{1},'SliceThickness'),
z = hdr{1}.SliceThickness;
else
z = 1;
end;
end;
vox = [hdr{1}.PixelSpacing z];
pos = hdr{1}.ImagePositionPatient';
dicom_to_patient = [orient*diag(vox) pos ; 0 0 0 1];
patient_to_tal = diag([-1 -1 1 1]);
mat = patient_to_tal*dicom_to_patient*analyze_to_dicom;
if strcmp(hdr{1}.Modality,'CT') && numel(hdr)>1
shear = (hdr{1}.ImagePositionPatient-hdr{2}.ImagePositionPatient) ...
* reshape(hdr{1}.ImageOrientationPatient,[3 2]);
if shear(1)
warning('shear(1) = %f not applied\n', shear(1));
end;
warning('shear(2) = %f applied\n', shear(2));
try
prms = spm_imatrix(mat);
prms(12) = shear(2);
mat = spm_matrix(prms);
end;
end;
% Possibly useful information
%-------------------------------------------------------------------
if checkfields(hdr{1},'AcquisitionTime','MagneticFieldStrength','MRAcquisitionType',...
'ScanningSequence','RepetitionTime','EchoTime','FlipAngle',...
'AcquisitionDate'),
tim = datevec(hdr{1}.AcquisitionTime/(24*60*60));
descrip = sprintf('%gT %s %s TR=%gms/TE=%gms/FA=%gdeg %s %d:%d:%.5g',...
hdr{1}.MagneticFieldStrength, hdr{1}.MRAcquisitionType,...
deblank(hdr{1}.ScanningSequence),...
hdr{1}.RepetitionTime,hdr{1}.EchoTime,hdr{1}.FlipAngle,...
datestr(hdr{1}.AcquisitionDate),tim(4),tim(5),tim(6));
else
descrip = hdr{1}.Modality;
end;
if ~true, % LEFT-HANDED STORAGE
mat = mat*[-1 0 0 (dim(1)+1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
end;
pinfo = [1 0 0]';
if isfield(hdr{1},'RescaleSlope') || isfield(hdr{1},'RescaleIntercept'),
pinfo = repmat(pinfo,1,length(hdr));
bytepix = spm_type('int16','bits')/8;
for i=1:length(hdr),
if isfield(hdr{i},'RescaleSlope'), pinfo(1,i) = hdr{i}.RescaleSlope; end;
if isfield(hdr{i},'RescaleIntercept'), pinfo(2,i) = hdr{i}.RescaleIntercept; end;
pinfo(3,i) = dim(1)*dim(2)*(i-1)*bytepix;
end;
end;
% Write the image volume
%-------------------------------------------------------------------
spm_progress_bar('Init',length(hdr),['Writing ' fname], 'Planes written');
%V = struct('fname',fname, 'dim',dim, 'dt',dt, 'pinfo',pinfo, 'mat',mat, 'descrip',descrip);
%V = spm_create_vol(V);
N = nifti;
N.dat = file_array(fname,dim,dt,0,pinfo(1),pinfo(2));
N.mat = mat;
N.mat0 = mat;
N.mat_intent = 'Scanner';
N.mat0_intent = 'Scanner';
N.descrip = descrip;
create(N);
for i=1:length(hdr),
plane = read_image_data(hdr{i});
if isfield(hdr{i},'RescaleSlope'), plane = plane*hdr{i}.RescaleSlope; end;
if isfield(hdr{i},'RescaleIntercept'), plane = plane+hdr{i}.RescaleIntercept; end;
plane = fliplr(plane);
if ~true, plane = flipud(plane); end; % LEFT-HANDED STORAGE
%V = spm_write_plane(V,plane,i);
N.dat(:,:,i) = plane;
spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
return;
%_______________________________________________________________________
%_______________________________________________________________________
function [images,guff] = select_tomographic_images(hdr)
images = {};
guff = {};
for i=1:length(hdr),
if ~checkfields(hdr{i},'Modality') || ~(strcmp(hdr{i}.Modality,'MR') ||...
strcmp(hdr{i}.Modality,'PT') || strcmp(hdr{i}.Modality,'CT'))
disp(['Cant find appropriate modality information for "' hdr{i}.Filename '".']);
guff = {guff{:},hdr{i}};
elseif ~checkfields(hdr{i},'StartOfPixelData','SamplesperPixel',...
'Rows','Columns','BitsAllocated','BitsStored','HighBit','PixelRepresentation'),
disp(['Cant find "Image Pixel" information for "' hdr{i}.Filename '".']);
guff = {guff{:},hdr{i}};
elseif ~checkfields(hdr{i},'PixelSpacing','ImagePositionPatient','ImageOrientationPatient'),
disp(['Cant find "Image Plane" information for "' hdr{i}.Filename '".']);
guff = {guff{:},hdr{i}};
elseif ~checkfields(hdr{i},'PatientID','InstanceNumber'),
%elseif ~checkfields(hdr{i},'PatientID','SeriesNumber','AcquisitionNumber','InstanceNumber'),
disp(['Cant find suitable filename info for "' hdr{i}.Filename '".']);
guff = {guff{:},hdr{i}};
else
images = {images{:},hdr{i}};
end;
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function [mosaic,standard] = select_mosaic_images(hdr)
mosaic = {};
standard = {};
for i=1:length(hdr),
if ~checkfields(hdr{i},'ImageType','CSAImageHeaderInfo') ||...
isfield(hdr{i}.CSAImageHeaderInfo,'junk') ||...
isempty(read_AcquisitionMatrixText(hdr{i})) ||...
isempty(read_NumberOfImagesInMosaic(hdr{i}))
standard = {standard{:},hdr{i}};
else
mosaic = {mosaic{:},hdr{i}};
end;
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function ok = checkfields(hdr,varargin)
ok = 1;
for i=1:(nargin-1),
if ~isfield(hdr,varargin{i}),
ok = 0;
break;
end;
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function clean = strip_unwanted(dirty)
msk = find((dirty>='a'&dirty<='z') | (dirty>='A'&dirty<='Z') |...
(dirty>='0'&dirty<='9') | dirty=='_');
clean = dirty(msk);
return;
%_______________________________________________________________________
%_______________________________________________________________________
function img = read_image_data(hdr)
img = [];
if hdr.SamplesperPixel ~= 1,
warning([hdr.Filename ': SamplesperPixel = ' num2str(hdr.SamplesperPixel) ' - cant be an MRI']);
return;
end;
prec = ['ubit' num2str(hdr.BitsAllocated) '=>' 'uint32'];
if isfield(hdr,'TransferSyntaxUID') && strcmp(hdr.TransferSyntaxUID,'1.2.840.10008.1.2.2') && strcmp(hdr.VROfPixelData,'OW'),
fp = fopen(hdr.Filename,'r','ieee-be');
else
fp = fopen(hdr.Filename,'r','ieee-le');
end;
if fp==-1,
warning([hdr.Filename ': cant open file']);
return;
end;
fseek(fp,hdr.StartOfPixelData,'bof');
img = fread(fp,hdr.Rows*hdr.Columns,prec);
fclose(fp);
if length(img)~=hdr.Rows*hdr.Columns,
error([hdr.Filename ': cant read whole image']);
end;
img = bitshift(img,hdr.BitsStored-hdr.HighBit-1);
if hdr.PixelRepresentation,
% Signed data - done this way because bitshift only
% works with signed data. Negative values are stored
% as 2s complement.
neg = logical(bitshift(bitand(img,uint32(2^hdr.HighBit)),-hdr.HighBit));
msk = (2^hdr.HighBit - 1);
img = double(bitand(img,msk));
img(neg) = img(neg)-2^(hdr.HighBit);
else
% Unsigned data
msk = (2^(hdr.HighBit+1) - 1);
img = double(bitand(img,msk));
end;
img = reshape(img,hdr.Columns,hdr.Rows);
return;
%_______________________________________________________________________
%_______________________________________________________________________
function nrm = read_SliceNormalVector(hdr)
str = hdr.CSAImageHeaderInfo;
val = get_numaris4_val(str,'SliceNormalVector');
for i=1:3,
nrm(i,1) = sscanf(val(i,:),'%g');
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function n = read_NumberOfImagesInMosaic(hdr)
str = hdr.CSAImageHeaderInfo;
val = get_numaris4_val(str,'NumberOfImagesInMosaic');
n = sscanf(val','%d');
if isempty(n), n=[]; end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function dim = read_AcquisitionMatrixText(hdr)
str = hdr.CSAImageHeaderInfo;
val = get_numaris4_val(str,'AcquisitionMatrixText');
dim = sscanf(val','%d*%d')';
if length(dim)==1,
dim = sscanf(val','%dp*%d')';
end;
if isempty(dim), dim=[]; end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function val = get_numaris4_val(str,name)
name = deblank(name);
val = {};
for i=1:length(str),
if strcmp(deblank(str(i).name),name),
for j=1:str(i).nitems,
if str(i).item(j).xx(1),
val = {val{:} str(i).item(j).val};
end;
end;
break;
end;
end;
val = strvcat(val{:});
return;
%_______________________________________________________________________
%_______________________________________________________________________
function fname = getfilelocation(hdr,root_dir)
switch root_dir
case 'projects_root'
dname = sprintf('/projects/%s%s%s%s%s_%.4d',...
strip_unwanted(strrep(hdr.PatientsName,'Medphys','')),...
filesep,...
strip_unwanted(hdr.StudyID), filesep,...
strip_unwanted(hdr.ProtocolName), ...
hdr.SeriesNumber);
case 'projects_cur'
dname = sprintf('%s%s%s%s%s%s%s_%.4d',...
pwd, filesep, ...
strip_unwanted(strrep(hdr.PatientsName,'Medphys','')),...
filesep,...
strip_unwanted(hdr.StudyID), filesep,...
strip_unwanted(hdr.ProtocolName), ...
hdr.SeriesNumber);
case 'patid',
dname = sprintf('%s%s%s%s%s_%.4d', pwd, filesep, ...
strip_unwanted(hdr.PatientID), ...
filesep, strip_unwanted(hdr.ProtocolName), ...
hdr.SeriesNumber);
case 'patname',
dname = sprintf('%s%s%s%s%s%s%s_%.4d', pwd, filesep, ...
strip_unwanted(hdr.PatientsName), ...
filesep, strip_unwanted(hdr.PatientID), ...
filesep, strip_unwanted(hdr.ProtocolName), ...
hdr.SeriesNumber);
otherwise
error('unknown file root specification');
end;
if exist(dname) ~= 7
mkdir_rec(dname);
end;
switch root_dir
case {'projects_root', 'projects_cur'},
fname = sprintf('f%s-%.5d-%.5d-%d.img', strip_unwanted(hdr.StudyID), ...
hdr.AcquisitionNumber,hdr.InstanceNumber, ...
hdr.EchoNumbers);
case {'patid', 'patname'},
fname = sprintf('f%s-%.5d-%.5d-%d.img', strip_unwanted(hdr.PatientID), ...
hdr.AcquisitionNumber,hdr.InstanceNumber, ...
hdr.EchoNumbers);
end;
fname = fullfile(dname, fname);
%_______________________________________________________________________
%_______________________________________________________________________
function suc = mkdir_rec(str);
% works on full pathnames only
if str(end) ~= filesep, str = [str filesep];end;
pos = findstr(str,filesep);
suc = zeros(1,length(pos));
for g=2:length(pos)
if exist(str(1:pos(g)-1)) ~= 7
suc(g) = mkdir(str(1:pos(g)-1));
end;
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function set_userdata(V,hdr)
% Get Diffusion information
try
userdata=struct('b', [], 'g', [], 'mat',[]);
userdata.b = str2num(get_numaris4_val(hdr.CSAImageHeaderInfo,'B_value'));
tmpg = str2num(strvcat(cellstr(get_numaris4_val(hdr.CSAImageHeaderInfo, ...
'DiffusionGradientDirection'))))';
if isempty(tmpg) && (userdata.b==0)
tmpg=[0 0 0];
end;
userdata.g = [tmpg(1) -tmpg(2) tmpg(3)];
catch
try
userdata=struct('b', [], 'g', [], 'mat',[]);
if strncmp(hdr.ImageComments,'No DW, scan',11)
userdata.b = 0;
userdata.g = [0 0 0];
elseif strncmp(hdr.ImageComments, 'DW: ', 4)
userdata.b = 1000; % fixed, not saved in images
[userdata.g(1) userdata.g(2) userdata.g(3)] = strread(hdr.ImageComments, ...
'DW: %f %f %f');
userdata.g(2) = -userdata.g(2); % this is empirical knowledge
% taken from above
else
userdata=[];
end
catch
userdata = [];
end;
end;
if ~isempty(userdata)
userdata.mat = V.mat;
spm_get_userdata(V.fname,userdata);
end;
%_______________________________________________________________________
%_______________________________________________________________________
function opts = spm_config_dicom
% Configuration file for dicom import jobs
%_______________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
% John Ashburner
% $Id: spm_config_dicom.m,v 1.4 2005/08/22 10:44:38 glauche Exp $
%_______________________________________________________________________
files.type = 'files';
files.name = 'DICOM files';
files.tag = 'files';
files.filter = 'any';
files.num = Inf;
files.help = {'Select the DICOM files to convert.'};
dirs.type = 'files';
dirs.name = 'DICOM directory tree(s)';
dirs.tag = 'dirs';
dirs.filter = 'dir';
dirs.num = [1 Inf];
dirs.help = {['Select the root(s) of the DICOM directory tree(s) to convert.'...
'The conversion will recursively search each directory for ' ...
'files to convert.']};
data.type = 'choice';
data.name = 'Select what (files or directories)?';
data.tag = 'data';
data.values = {files dirs};
data.help = {['Choose ''DICOM dir'' if your data is organised in a' ...
'directory tree. The conversion will then look recursively' ...
'in each subdirectory to collect the files']};
root.type = 'menu';
root.name = 'Root directory for converted images';
root.tag = 'root';
root.labels = {'/projects/<ProjectName>', './<ProjectName>', ...
'./<PatientID>', './<PatientName>'};
root.values = {'projects_root', 'projects_cur', ...
'patid', 'patname'};
root.def = 'util.dicom.root';
root.help = {'Choose root directory of converted file tree.'};
opts.type = 'branch';
opts.name = 'DICOM Import';
opts.tag = 'dicom';
opts.val = {data,root};
opts.prog = @convert_dicom;
opts.help = {[...
'DICOM Conversion. Most scanners produce files in DICOM format. '...
'Thios routine attempts to convert DICOM files into SPM compatible '...
'image volumes, which are written into the current directory. '...
'Note that not all flavours '...
'of DICOM can be handled, as DICOM is a very complicated format, and '...
'some scanner manufacturers use their own fields, which are not '...
'in the official documentation at http://medical.nema.org/']};
% add defaults, if not already set
global defaults
if isfield(defaults,'util')
if isfield(defaults.util,'dicom')
return;
end;
end;
defaults.util.dicom.root = 'projects_root';
%------------------------------------------------------------------------
%------------------------------------------------------------------------
function convert_dicom(job)
if isfield(job.data,'files')
files = {strvcat(job.data.files)};
else
for k = 1:numel(job.data.dirs)
files1=recurse_dirs(job.data.dirs{k},[]);
if k==1
files = files1;
else
files = {files1{:} files{:}};
end;
end;
end;
for k=1:numel(files)
try
hdr = spm_dicom_headers(strvcat(files{k}));
spm_dicom_convert(hdr,'all',job.root);
catch
[p n e v] = fileparts(deblank(files{k}(1,:)));
fprintf(['\n--------------------------\nerror '...
'converting images in directory\n\n%s.\n\n'...
' Error message was:\n%s\nSkipping to next directory\n'], p, lasterr);
hdr = [];
end;
end;
return;
%______________________________________________________________________
% Recursively get all files from a directory and all its subdirectories
function pp = recurse_dirs(cdir,pp)
pp1 = [];
ad = dir(cdir);
for i=1:length(ad)
if ad(i).isdir
if ~strcmp(ad(i).name,'.') & ~strcmp(ad(i).name,'..')
pp = recurse_dirs([cdir filesep ad(i).name],pp);
end;
else
pp1 = strvcat(pp1,[cdir filesep ad(i).name]);
end;
end;
pp{end+1} = pp1;
|
