Skip to content

Commit

Permalink
Second fix e0404#503
Browse files Browse the repository at this point in the history 
Made it so the interpolation is essential when a grid is present.
  • Loading branch information
@NathanKunz
NathanKunz committed May 5, 2021
1 parent 6ce1e78 commit 677dd90
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 21 deletions.
1 change: 0 additions & 1 deletion dicom/matRad_importDicomCt.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
% call
% ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool)
% ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid)
% ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, visBool)
% ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid, visBool)
%
% input
Expand Down
55 changes: 39 additions & 16 deletions dicom/matRad_interpDicomCtCube.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@

coordsOfFirstPixel = [origCtInfo.ImagePositionPatient];

% set up grid vectors
% set up original grid vectors
x = coordsOfFirstPixel(1,1) + origCtInfo(1).ImageOrientationPatient(1) * ...
origCtInfo(1).PixelSpacing(1)*double([0:origCtInfo(1).Columns-1]);
y = coordsOfFirstPixel(2,1) + origCtInfo(1).ImageOrientationPatient(5) * ...
Expand All @@ -59,29 +59,52 @@
yq(yq > max(yqRe)) = [];
zq(zq < min(zqRe)) = [];
zq(zq > max(zqRe)) = [];
else
xq = coordsOfFirstPixel(1,1):origCtInfo(1).ImageOrientationPatient(1)*resolution.x: ...
(coordsOfFirstPixel(1,1)+origCtInfo(1).ImageOrientationPatient(1)*origCtInfo(1).PixelSpacing(1)*double(origCtInfo(1).Columns-1));
yq = [coordsOfFirstPixel(2,1):origCtInfo(1).ImageOrientationPatient(5)*resolution.y: ...
(coordsOfFirstPixel(2,1)+origCtInfo(1).ImageOrientationPatient(5)*origCtInfo(1).PixelSpacing(2)*double(origCtInfo(1).Rows-1))];
zq = coordsOfFirstPixel(3,1):resolution.z: coordsOfFirstPixel(3,end);
end


% skip interpolation if the given resolution is equal to the orignal resolution
tmpCmp = abs([origCtInfo(1).PixelSpacing(1),origCtInfo(1).PixelSpacing(2),origCtInfo(1).SliceThickness] - [resolution.x,resolution.y,resolution.z]);
if(all(tmpCmp < 1e-6))
interpCt.cubeIV{1} = origCt;
else

% set up grid matrices - implicit dimension permuation (X Y Z-> Y X Z)
% Matlab represents internally in the first matrix dimension the
% ordinate axis and in the second matrix dimension the abscissas axis
[ X, Y, Z] = meshgrid(x,y,z);
[Xq, Yq, Zq] = meshgrid(xq,yq,zq);
% interpolate cube - cube is now stored in Y X Z
interpCt.cubeIV{1} = interp3(X,Y,Z,double(origCt),Xq,Yq,Zq);

else
% skip interpolation if the given resolution is equal (with an acceptance of 1e^-6) to the orignal resolution
if(all(abs([origCtInfo(1).PixelSpacing(1),origCtInfo(1).PixelSpacing(2),origCtInfo(1).SliceThickness] - [resolution.x,resolution.y,resolution.z]) < 1e-6))

% set the grid and resolution to the original values, to avoid reounding errors
resolution.x = origCtInfo(1).PixelSpacing(1);
resolution.y = origCtInfo(1).PixelSpacing(2);
resolution.z = origCtInfo(1).SliceThickness;
xq = x;
yq = y;
zq = z;

% set the return ctCube to the original Ct information because
% there is no interpolation needed
interpCt.cubeIV{1} = origCt;
cfg = MatRad_Config.instance();

else
% calculate new grid for the interpolation,
% grid equals a range of the first pixel, to the original
% last pixel, with steps equal the given resolution
xq = coordsOfFirstPixel(1,1):origCtInfo(1).ImageOrientationPatient(1)*resolution.x: ...
(coordsOfFirstPixel(1,1)+origCtInfo(1).ImageOrientationPatient(1)*origCtInfo(1).PixelSpacing(1)*double(origCtInfo(1).Columns-1));
yq = [coordsOfFirstPixel(2,1):origCtInfo(1).ImageOrientationPatient(5)*resolution.y: ...
(coordsOfFirstPixel(2,1)+origCtInfo(1).ImageOrientationPatient(5)*origCtInfo(1).PixelSpacing(2)*double(origCtInfo(1).Rows-1))];
zq = coordsOfFirstPixel(3,1):resolution.z: coordsOfFirstPixel(3,end);
% set up grid matrices - implicit dimension permuation (X Y Z-> Y X Z)
% Matlab represents internally in the first matrix dimension the
% ordinate axis and in the second matrix dimension the abscissas axis
[ X, Y, Z] = meshgrid(x,y,z);
[Xq, Yq, Zq] = meshgrid(xq,yq,zq);
% interpolate cube - cube is now stored in Y X Z
interpCt.cubeIV{1} = interp3(X,Y,Z,double(origCt),Xq,Yq,Zq);

end
end

% some meta information
interpCt.resolution = resolution;

Expand Down
10 changes: 6 additions & 4 deletions dicom/matRad_scanDicomImportFolder.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,8 @@
%% check for dicom files and differentiate patients, types, and series
numOfFiles = numel(fileList(:,1));
h = waitbar(0,'Please wait...');
% precision value for double to string conversion
str2numPrc = 10;
%h.WindowStyle = 'Modal';
steps = numOfFiles;
for i = numOfFiles:-1:1
Expand Down Expand Up @@ -95,7 +97,7 @@

try
if strcmp(info.Modality,'CT')
fileList{i,9} = num2str(info.PixelSpacing(1),15);
fileList{i,9} = num2str(info.PixelSpacing(1),str2numPrc);
else
fileList{i,9} = NaN;
end
Expand All @@ -104,7 +106,7 @@
end
try
if strcmp(info.Modality,'CT')
fileList{i,10} = num2str(info.PixelSpacing(2),15);
fileList{i,10} = num2str(info.PixelSpacing(2),str2numPrc);
else
fileList{i,10} = NaN;
end
Expand All @@ -116,9 +118,9 @@
%usually the Attribute should be SliceThickness, but it
%seems like some data uses "SpacingBetweenSlices" instead.
if isfield(info,'SliceThickness') && ~isempty(info.SliceThickness)
fileList{i,11} = num2str(info.SliceThickness,15);
fileList{i,11} = num2str(info.SliceThickness,str2numPrc);
elseif isfield(info,'SpacingBetweenSlices')
fileList{i,11} = num2str(info.SpacingBetweenSlices,15);
fileList{i,11} = num2str(info.SpacingBetweenSlices,str2numPrc);
else
matRad_cfg.dispError('Could not identify spacing between slices since neither ''SliceThickness'' nor ''SpacingBetweenSlices'' are specified');
end
Expand Down

0 comments on commit 677dd90

Please sign in to comment.