Hi Darren, Mateus & list,
you might want to try the attached files instead of the SPM5 original
ones. There is an option for dicom conversion into .nii or .hdr/.img.
Please let me know whether they work, I have just assembled them from our
in-house version without additional testing.
Volkmar
On Mon, 4 Dec 2006, Darren G wrote:
> you will now have a nifti file and the hdr+img files can be deleted. You
could put this into a loop
> and convert all the converted dicom files to nii, backup and delete the
original hdr+img files and
> then according to what Volkmar said must of the rest of the output would
be .nii files.
>
> Darren
>
> ==============Original message text===============
> On Mon, 04 Dec 2006 4:39:19 pm CST Volkmar Glauche wrote:
>
> Actually, for most parts in SPM output is .nii if your input is in .nii
> format.
> Personally, I favour the single file format, because it allows to optimise the
> file system for larger files without wasting a lot of disk space for just a
> small header file.
>
> Zitat von "Ashburner John (PSYCHOLOGY)"
> <[log in to unmask]>:
>
> > There's no easy way of doing this without heavily modifying SPM.
> >
> > The reason why SPM5 produces .img and .hdr files instead of .nii files
> > is that we have a backup system at the FIL. Files that are modified are
> > backed up more frequently, and I figured it would not be such a good
> > idea to have the entire .nii file backed up whenever a little
> > modification is made to a header.
> >
> > Best regards,
> > -John
> >
> > -----Original Message-----
> > From: SPM (Statistical Parametric Mapping) [mailto:[log in to unmask]]
> > On Behalf Of Mateus Joffily
> > Sent: Monday, December 04, 2006 8:34 PM
> > To: [log in to unmask]
> > Subject: [SPM] .nii files
> >
> > Hi,
> >
> > Is there a way of forcing spm5 to write the image files with '.nii'
> > extension, instead of '.hdr+.img'?
> >
> > Thanks for your help,
> > Mateus
> >
> >
> ===========End of original message text===========
>
>
>
> --------------------------------
> Darren R. Gitelman, M.D.
> Cognitive Neurology and the Alzheimer's Disease Center
> Northwestern Univ., 320 E. Superior St., Searle 11-470, Chicago, IL 60611
> Voice: (312) 908-9023 Fax: (312) 908-8789
> --------------------------------
>
>
--
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 http://fbi.uniklinik-freiburg.de/
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 531 2006-05-16 09:17:55Z volkmar $
%_______________________________________________________________________
data.type = 'files';
data.name = 'DICOM files';
data.tag = 'data';
data.filter = '.*';
data.num = Inf;
data.help = {'Select the DICOM files to convert.'};
outdir.type = 'files';
outdir.name = 'Output directory';
outdir.tag = 'outdir';
outdir.filter = 'dir';
outdir.num = 1;
outdir.val = {''};
outdir.help = {[...
'Select a directory where files are written. '...
'Default is current directory.']};
root.type = 'menu';
root.name = 'Directory structure for converted files';
root.tag = 'root';
root.labels = {'Output directory: ./<StudyDate-StudyTime>', ...
'Output directory: ./<PatientID>', ...
'Output directory: ./<PatientID>/<StudyDate-StudyTime>', ...
'Output directory: ./<PatientName>', ...
'No directory hierarchy'};
root.values = {'date_time',...
'patid', 'patid_date', 'patname','flat'};
root.def = 'util.dicom.root';
root.help = {['Choose root directory of converted file tree. The options ' ...
'are:'], '',...
['* Output directory: ./<StudyDate-StudyTime>: ' ...
'Automatically determine the project name and try to ' ...
'convert into the output directory, starting with '...
'a StudyDate-StudyTime subdirectory. This option is useful if automatic '...
'project recognition fails and one wants to convert data into '...
'a project directory.'], '',...
['* Output directory: ./<PatientID>: ' ...
'Convert into the output directory, starting with '...
'a PatientID subdirectory.'], '',...
['* Output directory: ./<PatientName>: ' ...
'Convert into the output directory, starting with '...
'a PatientName subdirectory.'],...
['* No directory hierarchy: Convert all files into the output ' ...
'directory, without sequence/series subdirectories']};
format.type = 'menu';
format.name = 'Output image format';
format.tag = 'format';
format.labels = {'Two file (img+hdr) NIfTI', 'Single file (nii) NIfTI'};
format.values = {'img', 'nii'};
format.val = {'nii'};
format.help = {['DICOM conversion can create separate img and hdr files ' ...
'or combine them in one file. The single file option will ' ...
'help you save space on your harddisk, but may be ' ...
'incompatible with programs that are not NIfTI-aware.'],...
['In any case, only 3D image files will be produced.']};
icedims.type = 'menu';
icedims.name = 'Use ICEDims in filename';
icedims.tag = 'icedims';
icedims.labels = {'No','Yes'};
icedims.values = {0, 1};
icedims.val = {0};
icedims.help = {['If image sorting fails, one can try using the additional ' ...
'SIEMENS ICEDims information to create unique filenames. ' ...
'Use this only if there would be multiple volumes with '...
'exactly the same file names.']};
convopts.type = 'branch';
convopts.name = 'Conversion options';
convopts.tag = 'convopts';
convopts.val = {format,icedims};
opts.type = 'branch';
opts.name = 'DICOM Import';
opts.tag = 'dicom';
opts.val = {data,root,outdir,convopts};
opts.prog = @convert_dicom;
opts.help = {[...
'DICOM Conversion. Most scanners produce data in DICOM format. '...
'This routine attempts to convert DICOM files into SPM compatible '...
'image volumes, which are written into the current directory by '...
'default. 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 = 'flat';
%------------------------------------------------------------------------
%------------------------------------------------------------------------
function convert_dicom(job)
wd = pwd;
try
if ~isempty(job.outdir{:})
cd(job.outdir{:});
fprintf(' Changing directory to: %s\n', job.outdir{:});
end
catch
error('Failed to change directory. Aborting DICOM import.');
end
if job.convopts.icedims
root_dir = ['ice' job.root];
else
root_dir = job.root;
end;
hdr = spm_dicom_headers(strvcat(job.data));
spm_dicom_convert(hdr,'all',root_dir,job.convopts.format);
if ~isempty(job.outdir)
fprintf(' Changing back to directory: %s\n', wd);
cd(wd);
end
function spm_dicom_convert(hdr,opts,root_dir,format)
% Convert DICOM images into something that SPM can use
% FORMAT spm_dicom_convert(hdr,opts,root_dir,format)
% 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
% 'spect' - SIEMENS Spectroscopy DICOMs (position only)
% This will write out a mask image volume with 1's
% set at the position of spectroscopy voxel(s).
% 'raw' - convert raw FIDs (not implemented)
% root_dir - 'flat' - SPM5 standard, do not produce file tree
% With all other options, files will be sorted into
% directories according to their sequence/protocol names
% 'date_time' - Place files under ./<StudyDate-StudyTime>
% 'patid' - Place files under ./<PatID>
% 'patid_date' - Place files under ./<PatID-StudyDate>
% 'name' - Place files under ./<PatName>
% format - output format
% 'img' Two file (hdr+img) NIfTI format
% 'nii' Single file NIfTI format
% All images will contain a single 3D dataset, 4D images
% will not be created.
%_______________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
% John Ashburner & Jesper Andersson
% $Id: spm_dicom_convert.m 698 2006-11-24 20:12:08Z volkmar $
if nargin<2, opts = 'all'; end;
if nargin<3, root_dir='flat';end;
if nargin<4, format='img';end;
[images,guff] = select_tomographic_images(hdr);
[spect,images] = select_spectroscopy_images(images);
[mosaic,standard] = select_mosaic_images(images);
if (strcmp(opts,'all') || strcmp(opts,'mosaic')) && ~isempty(mosaic),
convert_mosaic(mosaic,root_dir,format);
end;
if (strcmp(opts,'all') || strcmp(opts,'standard')) && ~isempty(standard),
convert_standard(standard,root_dir,format);
end;
if (strcmp(opts,'all') || strcmp(opts,'spect')) && ~isempty(spect),
convert_spectroscopy(spect,root_dir,format);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function convert_mosaic(hdr,root_dir,format)
spm_progress_bar('Init',length(hdr),'Writing Mosaic', 'Files written');
for i=1:length(hdr),
% Output filename
%-------------------------------------------------------------------
fname = getfilelocation(hdr{i},root_dir,'f',format);
% 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 = N.dat;
dat.scl_slope = [];
dat.scl_inter = [];
for j=1:dim(3),
dat(:,:,j) = volume(:,:,j);
end;
spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
return;
%_______________________________________________________________________
%_______________________________________________________________________
function convert_standard(hdr,root_dir,format)
hdr = sort_into_volumes(hdr);
for i=1:length(hdr),
write_volume(hdr{i},root_dir,format);
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);
% This line is a fudge because of some problematic data that Bogdan,
% Cynthia and Stefan were trying to convert. I hope it won't cause
% problems for others -JA
dist2 = 0;
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 there 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 &&...
sum((hdr{i}.ImageOrientationPatient - vol{j}{1}.ImageOrientationPatient).^2)<1e-5 &&...
sum((hdr{i}.PixelSpacing - vol{j}{1}.PixelSpacing).^2)<1e-5 && ...
identical_ice_dims && dist2<1e-3;
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;
% For raw image data, tell apart real/complex or phase/magnitude
if isfield(hdr{i},'ImageType') && isfield(vol{j}{1}, ...
'ImageType')
match = match && strcmp(hdr{i}.ImageType, ...
vol{j}{1}.ImageType);
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};
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)>1e-4,
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.4g *\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,format)
% Output filename
%-------------------------------------------------------------------
fname = getfilelocation(hdr{1}, root_dir,'s',format);
% 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 convert_spectroscopy(hdr,root_dir,format)
for i=1:length(hdr),
write_spectroscopy_volume(hdr(i),root_dir,format);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function write_spectroscopy_volume(hdr,root_dir,format)
% Output filename
%-------------------------------------------------------------------
fname = getfilelocation(hdr{1}, root_dir,'S',format);
% Image dimensions
%-------------------------------------------------------------------
nc = get_numaris4_numval(hdr{1}.Private_0029_1210,'Columns');
nr = get_numaris4_numval(hdr{1}.Private_0029_1210,'Rows');
dim = [nc nr numel(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(get_numaris4_numval(hdr{1}.Private_0029_1210,...
'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) = get_numaris4_numval(hdr{i}.Private_0029_1210,...
'ImagePositionPatient')*orient(:,3);
end;
z = mean(diff(z));
else
try
z = get_numaris4_numval(hdr{1}.Private_0029_1210,...
'SliceThickness');
catch
z = 1;
end;
end;
vox = [get_numaris4_numval(hdr{1}.Private_0029_1210,'PixelSpacing') z];
pos = get_numaris4_numval(hdr{1}.Private_0029_1210,'ImagePositionPatient')';
%dicom_to_patient = [orient*diag(vox) pos-1.5*orient*([0 vox(2) 0]') ; 0 0 0 1];
dicom_to_patient = [orient*diag(vox) pos ; 0 0 0 1];
patient_to_tal = diag([-1 -1 1 1]);
warning('Don''t know exactly what positions in spectroscopy files should be - just guessing!')
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_spect_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')||isfield(hdr{i},'Private_7fe1_0010')),
disp(['Cant find "Image Pixel" information for "' hdr{i}.Filename '".']);
guff = {guff{:},hdr{i}};
elseif ~(checkfields(hdr{i},'PixelSpacing','ImagePositionPatient','ImageOrientationPatient')||isfield(hdr{i},'Private_0029_1210')),
disp(['Cant find "Image Plane" information for "' hdr{i}.Filename '".']);
guff = {guff{:},hdr{i}};
elseif ~checkfields(hdr{i},'PatientID','SeriesNumber','AcquisitionNumber','InstanceNumber'),
disp(['Cant find suitable filename info for "' hdr{i}.Filename '".']);
if ~isfield(hdr{i},'SeriesNumber')
disp('Setting SeriesNumber to 1');
hdr{i}.SeriesNumber=1;
images = {images{:},hdr{i}};
end;
if ~isfield(hdr{i},'AcquisitionNumber')
disp('Setting AcquisitionNumber to 1');
hdr{i}.AcquisitionNumber=1;
images = {images{:},hdr{i}};
end;
if ~isfield(hdr{i},'InstanceNumber')
disp('Setting InstanceNumber to 1');
hdr{i}.InstanceNumber=1;
images = {images{:},hdr{i}};
end;
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 [spect,images] = select_spectroscopy_images(hdr)
spectsel = zeros(1,numel(hdr));
for i=1:length(hdr),
if isfield(hdr{i},'SOPClassUID')
spectsel(i) = strcmp(hdr{i}.SOPClassUID,'1.3.12.2.1107.5.9.1');
end;
end;
spect = hdr(logical(spectsel));
images= hdr(~logical(spectsel));
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 img = read_spect_data(hdr)
% Image dimensions
%-------------------------------------------------------------------
nc = get_numaris4_numval(hdr.Private_0029_1210,'Columns');
nr = get_numaris4_numval(hdr.Private_0029_1210,'Rows');
img = ones(nr,nc);
%_______________________________________________________________________
%_______________________________________________________________________
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 val = get_numaris4_numval(str,name)
val1 = get_numaris4_val(str,name);
for k = 1:size(val1,1)
val(k)=str2num(val1(k,:));
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function fname = getfilelocation(hdr,root_dir,prefix,format)
if nargin < 3
prefix = 'f';
end;
if strncmp(root_dir,'ice',3)
root_dir = root_dir(4:end);
imtype = textscan(hdr.ImageType,'%s','delimiter','\\');
try
imtype = imtype{1}{3};
catch
imtype = '';
end;
prefix = [prefix imtype get_numaris4_val(hdr.CSAImageHeaderInfo,'ICE_Dims')];
end;
if strcmp(root_dir, 'flat')
% Standard SPM5 file conversion
%-------------------------------------------------------------------
if checkfields(hdr,'SeriesNumber','AcquisitionNumber')
if checkfields(hdr,'EchoNumbers')
fname = sprintf('%s%s-%.4d-%.5d-%.6d-%.2d.%s', prefix, strip_unwanted(hdr.PatientID),...
hdr.SeriesNumber, hdr.AcquisitionNumber, hdr.InstanceNumber,...
hdr.EchoNumbers, format);
else
fname = sprintf('%s%s-%.4d-%.5d-%.6d.%s', prefix, strip_unwanted(hdr.PatientID),...
hdr.SeriesNumber, hdr.AcquisitionNumber, ...
hdr.InstanceNumber, format);
end;
else
fname = sprintf('%s%s-%.6d.%s',prefix, ...
strip_unwanted(hdr.PatientID),hdr.InstanceNumber, format);
end;
fname = fullfile(pwd,fname);
return;
end;
% more fancy stuff - sort images into subdirectories
if ~isfield(hdr,'ProtocolName')
if isfield(hdr,'SequenceName')
hdr.ProtocolName = hdr.SequenceName;
else
hdr.ProtocolName='unknown';
end;
end;
if ~isfield(hdr,'SeriesDescription')
hdr.SeriesDescription = 'unknown';
end;
if ~isfield(hdr,'EchoNumbers')
hdr.EchoNumbers = 0;
end;
m = sprintf('%02d', floor(rem(hdr.StudyTime/60,60)));
h = sprintf('%02d', floor(hdr.StudyTime/3600));
studydate = sprintf('%s_%s-%s', datestr(hdr.StudyDate,'yyyy-mm-dd'), ...
h,m);
serdes = strrep(strip_unwanted(hdr.SeriesDescription),...
strip_unwanted(hdr.ProtocolName),'');
protname = sprintf('%s%s_%.4d',strip_unwanted(hdr.ProtocolName), ...
serdes, hdr.SeriesNumber);
switch root_dir
case 'date_time',
root = pwd;
dname = fullfile(root, studydate, protname);
case 'patid',
dname = fullfile(pwd, strip_unwanted(hdr.PatientID), ...
protname);
case 'patid_date',
dname = fullfile(pwd, strip_unwanted(hdr.PatientID), ...
studydate, protname);
case 'patname',
dname = fullfile(pwd, strip_unwanted(hdr.PatientsName), ...
strip_unwanted(hdr.PatientID), ...
protname);
otherwise
error('unknown file root specification');
end;
if ~exist(dname,'dir'),
mkdir_rec(dname);
end;
switch root_dir
case 'date_time',
fname = sprintf('%s%s-%.5d-%.5d-%d.%s', prefix, studydate, ...
hdr.AcquisitionNumber,hdr.InstanceNumber, ...
hdr.EchoNumbers,format);
case {'patid', 'patid_date', 'patname'},
fname = sprintf('%s%s-%.5d-%.5d-%d.%s', prefix, strip_unwanted(hdr.PatientID), ...
hdr.AcquisitionNumber,hdr.InstanceNumber, ...
hdr.EchoNumbers,format);
end;
fname = fullfile(dname, fname);
%_______________________________________________________________________
%_______________________________________________________________________
function suc = mkdir_rec(str)
% works on full pathnames only
opwd=pwd;
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),'dir'),
cd(str(1:pos(g-1)-1));
suc(g) = mkdir(str(pos(g-1)+1:pos(g)-1));
end;
end;
cd(opwd);
return;
%_______________________________________________________________________
%_______________________________________________________________________
function ret = read_ascconv(hdr)
% In SIEMENS data, there is an ASCII text section with
% additional information items. This section starts with a code
% ### ASCCONV BEGIN ###
% and ends with
% ### ASCCONV END ###
% It is read by spm_dicom_headers into an entry 'MrProtocol' in CSASeriesHeaderInfo
% The additional items are assignments in C syntax, here they are just
% translated according to
% [] -> ()
% " -> '
% 0xX -> hex2dec('X')
% and collected in a struct.
ret=struct;
% get ascconv data
X=get_numaris4_val(hdr.CSASeriesHeaderInfo,'MrProtocol');
ascstart = findstr(X,'### ASCCONV BEGIN ###');
ascend = findstr(X,'### ASCCONV END ###');
if ~isempty(ascstart) && ~isempty(ascend)
tokens = textscan(char(X((ascstart+22):(ascend-1))),'%s', ...
'delimiter',char(10));
for k = 1:numel(tokens{1})
tokens{1}{k}=regexprep(tokens{1}{k},'\[([0-9]*)\]','($1+1)');
tokens{1}{k}=regexprep(tokens{1}{k},'"(.*)"','''$1''');
tokens{1}{k}=regexprep(tokens{1}{k},'0x([0-9a-fA-F]*)','hex2dec(''$1'')');
try
eval(['ret.' tokens{1}{k} ';']);
catch
disp(['AscConv: Error evaluating ''ret.' tokens{1}{k} ''';']);
end;
end;
end;
%_______________________________________________________________________
|
