Homec4science
Diffusion segmentation (master) Browse

segmentation/functions/readDicomStack.m
01b7952df4bcmaster

Last Change
Raw File

readDicomStack.m
View Options

function [img, imgInfo] = readDicomStack(dicomDir, opt)
%
dicomdict('set','dicom.dic');
% kernelTypes = {'bone', 'lung', 'none'};
if(~isempty(dicomDir) && dicomDir(end) ~= filesep)
dicomDir = [dicomDir filesep];
end
if (nargin == 0)
dicomDir = './';
opt = 'max';
else
if (nargin == 1)
opt = 'max';
end
end
if class(dicomDir)=='char'
filesList = dir([dicomDir '*.dcm']);
if(isempty(filesList))
filesList = dir(dicomDir);
filesList = filesList(3:end);
end
tmp{size(filesList, 1)} = '';
for ii = 1:size(filesList, 1)
tmp{ii} = [dicomDir filesList(ii).name];
end
filesList = tmp;
else
filesList = dicomDir;
end
Headers{size(filesList, 2)} = [];
sliceLoc = zeros(size(filesList, 1), 1);
kernelType = zeros(size(sliceLoc));
for ii = 1:size(filesList, 2)
try
Headers{ii} = dicominfo([filesList{ii}]);
catch
sliceLoc(ii) = inf;
continue;
end
try
sliceLoc(ii) = Headers{ii}.ImagePositionPatient(3);
catch
sliceLoc(ii) = -inf;
end
if (isfield(Headers{ii}, 'ConvolutionKernel'))
switch lower(Headers{ii}.ConvolutionKernel)
case 'bone'
kernelType(ii) = 1;
case 'lung'
kernelType(ii) = 2;
case 'standard'
kernelType(ii) = 2;
case 'soft'
kernelType(ii) = 2;
otherwise
kernelType(ii) = 3;
end
else
kernelType(ii) = 3;
end
end
[kernelType, ind] = sort(kernelType);
sliceLoc = sliceLoc(ind);
Headers = Headers(:, ind);
numStacks = 0;
ii = 1;
while (ii<size(Headers, 2))
jj = ii;
while(jj<=size(Headers, 2) && kernelType(jj) == kernelType(ii)), jj = jj+1; end
tmp = sliceLoc(ii:jj-1);
[sliceLoc(ii:jj-1), ind] = sort(tmp);
tmp = Headers(:, ii:jj-1);
Headers(:, ii:jj-1) = tmp(:, ind);
ii = jj;
numStacks = numStacks+1;
end
if (numStacks == 0)
imgInfo = [];
img = [];
return;
end
currentSlice = 1;
img{numStacks} = [];
for ii = 1:numStacks
currentKernel = kernelType(currentSlice);
while(sliceLoc(currentSlice)==-inf), currentSlice = currentSlice+1; end
xMin = 0; xMax = 0; xSize = 0;
yMin = 0; yMax = 0; ySize = 0;
zMin = 0; zMax = 0; zSize = 0;
PixelSpacing = [0, 0];
SliceThickness = 0;
if (isfield(Headers{currentSlice}, 'PixelSpacing'))
PixelSpacing = Headers{currentSlice}.PixelSpacing;
end
% if (isfield(Headers{currentSlice}, 'SliceThickness'))
% SliceThickness = Headers{currentSlice}.SliceThickness;
% end
if (isfield(Headers{currentSlice}, 'ImagePositionPatient'))
xMin = Headers{currentSlice}.ImagePositionPatient(1);
yMin = Headers{currentSlice}.ImagePositionPatient(2);
zMin = Headers{currentSlice}.ImagePositionPatient(3);
end
if (isfield(Headers{currentSlice}, 'Columns'))
xSize = Headers{currentSlice}.Columns;
xMax = xMin + double(xSize - 1) * PixelSpacing(2);
end
if (isfield(Headers{currentSlice}, 'Rows'))
ySize = Headers{currentSlice}.Rows;
yMax = yMin + double(ySize - 1) * PixelSpacing(1);
end
switch currentKernel
case 1
imgInfo(ii).ConvolutionKernel = 'Bone';
case 2
imgInfo(ii).ConvolutionKernel = 'Lung';
otherwise
imgInfo(ii).ConvolutionKernel = 'None';
end
while(currentSlice <= numel(kernelType) && currentKernel == kernelType(currentSlice))
zSize = zSize+1;
try
tmp = dicomread(Headers{currentSlice});
if (isfield(Headers{currentSlice}, 'RescaleSlope') && Headers{currentSlice}.RescaleSlope~=0)
Slope = Headers{currentSlice}.RescaleSlope;
else
Slope = 1;
end
if (isfield(Headers{currentSlice}, 'RescaleIntercept'))
Intercept = Headers{currentSlice}.RescaleIntercept;
else
Intercept = 0;
end
catch
currentSlice = currentSlice + 1;
continue;
end
img{ii}(:, :, zSize) = tmp*Slope + Intercept;
if(SliceThickness == 0 && currentSlice<numel(kernelType) && kernelType(currentSlice+1) == currentKernel)
SliceThickness = sliceLoc(currentSlice+1) - sliceLoc(currentSlice);
else
newLoc = zMin+(zSize-1)*SliceThickness;
if(SliceThickness ~= 0 && newLoc ~= sliceLoc(currentSlice) && currentSlice<numel(kernelType) && kernelType(currentSlice+1) == currentKernel)
sliceLoc(currentSlice) = newLoc;
end
end
zMax = sliceLoc(currentSlice);
currentSlice = currentSlice + 1;
end
imgInfo(ii).BoundingBox.xMin = xMin; imgInfo(ii).BoundingBox.xMax = xMax; imgInfo(ii).xSize = xSize;
imgInfo(ii).BoundingBox.yMin = yMin; imgInfo(ii).BoundingBox.yMax = yMax; imgInfo(ii).ySize = ySize;
imgInfo(ii).BoundingBox.zMin = zMin; imgInfo(ii).BoundingBox.zMax = zMax; imgInfo(ii).zSize = zSize;
imgInfo(ii).SliceThickness = (zMax -zMin)/(zSize-1);
imgInfo(ii).PixelSpacing = [(xMax - xMin)/double(xSize-1) (yMax - yMin)/double(ySize-1)];
end
ind = [];
display(['TEST : length(img) = ' num2str(length(img))]);
display(['TEST : numStacks = ' num2str(numStacks)]);
display(['TEST : opt = ' num2str(opt)]);
for ii = 1:numStacks
if(isempty(img{ii}))
ind = [ind, ii];
end
end
img(:, ind) = [];
imgInfo(:, ind) = [];
switch opt
case 'all'
case 'first'
img = img{1};
imgInfo = imgInfo(1);
case 'last'
img = img{end};
imgInfo = imgInfo(end);
otherwise
SZs = zeros(size(imgInfo, 2), 1);
for tmp = 1:size(imgInfo, 2)
SZs(tmp) = imgInfo(tmp).zSize;
end
tmp = max(SZs);
ind = find(SZs == tmp);
if(numel(ind)>1)
for ii=1:numel(ind)
if(imgInfo(ind(ii)).ConvolutionKernel == 'Lung')
break;
end
end
if (ii>numel(ind))
ind = ind(1);
else
ind = ind(ii);
end
end
display(['TEST : length(img) = ' num2str(length(img))]);
display(['TEST : ind = ' num2str(ind)]);
imgInfo = imgInfo(ind);
end
% figure
% for ii = 1:size(img{1}, 3)
% imshow(uint8(img{1}(:, :, ii)/4095*255));
% pause(0.1)
% end
% % % stackNo = 1;
% % %
% % %
% % %
% % % info = dicominfo([dicomDir filesList(1).name]);
% % % tmp = dicomread(info);
% % % imgInfo(stackNo).BoundingBox.xMin = info.ImagePositionPatient(1);
% % % imgInfo(stackNo).BoundingBox.yMin = info.ImagePositionPatient(2);
% % % imgInfo(stackNo).BoundingBox.zMin = info.ImagePositionPatient(3);
% % % imgInfo(stackNo).PixelSpacing = info.PixelSpacing;
% % % imgInfo(stackNo).SliceThickness = 0;
% % % imgInfo(stackNo).xSize = info.Columns;
% % % imgInfo(stackNo).ySize = info.Rows;
% % % imgInfo(stackNo).zSize = 1;
% % % imgInfo(stackNo).BoundingBox.xMax = imgInfo(stackNo).BoundingBox.xMin + double(imgInfo(stackNo).xSize -1) * imgInfo(stackNo).PixelSpacing(1);
% % % imgInfo(stackNo).BoundingBox.yMax = imgInfo(stackNo).BoundingBox.yMin + double(imgInfo(stackNo).ySize -1) * imgInfo(stackNo).PixelSpacing(2);
% % % imgInfo(stackNo).BoundingBox.zMax = double(imgInfo(stackNo).BoundingBox.zMin);
% % % img{stackNo} = tmp;
% % % infoPre = info;
% % %
% % % numSlicesInStack = [];
% % % for ii = 2:size(filesList, 1)
% % % try
% % % info = dicominfo([dicomDir filesList(ii).name]);
% % % catch
% % % continue;
% % % end
% % % if ( ~isequal(info.Columns, infoPre.Columns) || ~isequal(info.Rows, infoPre.Rows) || ~isequal(info.PatientOrientation, infoPre.PatientOrientation))
% % % img{stackNo} = eval([class(tmp) '(zeros(imgInfo(stackNo).xSize, imgInfo(stackNo).ySize, imgInfo(stackNo).zSize))']);
% % % numSlicesInStack(stackNo) = imgInfo(stackNo).zSize;
% % % stackNo = stackNo + 1;
% % % try
% % % imgInfo(stackNo).BoundingBox.xMin = info.ImagePositionPatient(1);
% % % catch
% % % imgInfo(stackNo).BoundingBox.xMin = 0;
% % % end
% % % try
% % % imgInfo(stackNo).BoundingBox.yMin = info.ImagePositionPatient(2);
% % % catch
% % % imgInfo(stackNo).BoundingBox.yMin = 0;
% % % end
% % % try
% % % imgInfo(stackNo).BoundingBox.zMin = info.ImagePositionPatient(3);
% % % catch
% % % imgInfo(stackNo).BoundingBox.zMin = 0;
% % % end
% % % try
% % % imgInfo(stackNo).PixelSpacing = info.PixelSpacing;
% % % catch
% % % imgInfo(stackNo).PixelSpacing = 1;
% % % end
% % %
% % % imgInfo(stackNo).SliceThickness = 0;
% % % try
% % % imgInfo(stackNo).xSize = info.Columns;
% % % imgInfo(stackNo).ySize = info.Rows;
% % % catch
% % % tmp = dicomread(info);
% % % imgInfo(stackNo).xSize = size(tmp, 2);
% % % imgInfo(stackNo).ySize = size(tmp, 1);
% % % end
% % %
% % %
% % % imgInfo(stackNo).zSize = 1;
% % % imgInfo(stackNo).BoundingBox.xMax = imgInfo(stackNo).BoundingBox.xMin + (imgInfo(stackNo).xSize -1) * imgInfo(stackNo).PixelSpacing(1);
% % % imgInfo(stackNo).BoundingBox.yMax = imgInfo(stackNo).BoundingBox.yMin + (imgInfo(stackNo).ySize -1) * imgInfo(stackNo).PixelSpacing(2);
% % % imgInfo(stackNo).BoundingBox.zMax = imgInfo(stackNo).BoundingBox.zMin;
% % % infoPre = info;
% % % else
% % % try
% % % imgInfo(stackNo).BoundingBox.zMax = info.ImagePositionPatient(3);
% % % catch
% % % imgInfo(stackNo).BoundingBox.zMax = imgInfo(stackNo).BoundingBox.zMax + 1;
% % % end
% % % imgInfo(stackNo).zSize = imgInfo(stackNo).zSize + 1;
% % % imgInfo(stackNo).SliceThickness = (imgInfo(stackNo).BoundingBox.zMax - imgInfo(stackNo).BoundingBox.zMin)/imgInfo(stackNo).zSize;
% % % end
% % % end
% % % img{stackNo} = eval([class(tmp) '(zeros(imgInfo(stackNo).xSize, imgInfo(stackNo).ySize, imgInfo(stackNo).zSize))']);
% % % numSlicesInStack(stackNo) = imgInfo(stackNo).zSize;
% % %
% % % switch opt
% % % case 'all'
% % % f = 1;
% % % for ii = 1:length(numSlicesInStack)
% % % for jj = 1:numSlicesInStack(ii)
% % % info = dicominfo([dicomDir filesList(f).name]);
% % % tmp = dicomread(info);
% % % img{ii}(:, :, jj) = tmp;
% % % f = f + 1;
% % % end
% % % end
% % % case 'first'
% % % for ii = 1:numSlicesInStack(1)
% % % info = dicominfo([dicomDir filesList(ii).name]);
% % % tmp = dicomread(info);
% % % img{1}(:, :, ii) = tmp;
% % % end
% % % imgInfo = imgInfo(1);
% % % img = img{1};
% % % case 'last'
% % % if numel(numSlicesInStack) == 1
% % % fileInd = 0;
% % % else
% % % fileInd = sum(numSlicesInStack(end-1));
% % % end
% % % for ii = 1:numSlicesInStack(end)
% % % info = dicominfo([dicomDir filesList(ii+fileInd).name]);
% % % tmp = dicomread(info);
% % % img{end}(:, :, ii) = tmp;
% % % end
% % % imgInfo = imgInfo(end);
% % % img = img{end};
% % % otherwise
% % % [~, ind] = max(numSlicesInStack);
% % % if ind == 1
% % % fileInd = 0;
% % % else
% % % fileInd = sum(numSlicesInStack(1:ind-1));
% % % end
% % % for ii = 1:numSlicesInStack(ind)
% % % info = dicominfo([dicomDir filesList(ii+fileInd).name]);
% % % tmp = dicomread(info);
% % % img{ind}(:, :, ii) = tmp;
% % % end
% % % imgInfo = imgInfo(ind);
% % % img = img{ind};
% % % end
c4science · Help