Skip to content

Commit

Permalink
BUG: itkDCMTKFileReader Philips multi-frame MRI Z spacing
Browse files Browse the repository at this point in the history
Addresses #4131, Z spacing on an Philips MRI volume saved as MRExtendedStorage.
  • Loading branch information
hossbach-cgm authored and thewtex committed Mar 7, 2024
1 parent aa65f4b commit dbd46b4
Show file tree
Hide file tree
Showing 4 changed files with 236 additions and 33 deletions.
88 changes: 56 additions & 32 deletions Modules/IO/DCMTK/src/itkDCMTKFileReader.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -1172,7 +1172,12 @@ DCMTKFileReader::GetDimensions(unsigned short & rows, unsigned short & columns)
int
DCMTKFileReader::GetSpacing(double * const spacing) const
{
double _spacing[3];
spacing[0] = 1.0;
spacing[1] = 1.0;
spacing[2] = 1.0;

double _spacing[3] = { 1.0, 1.0, 1.0 };
bool _have_val[3] = { false, false, false };
//
// There are several tags that can have spacing, and we're going
// from most to least desirable, starting with PixelSpacing, which
Expand All @@ -1182,9 +1187,16 @@ DCMTKFileReader::GetSpacing(double * const spacing) const

// first, shared function groups sequence, then
// per-frame groups sequence
_spacing[0] = _spacing[1] = _spacing[2] = 0.0;
int rval(EXIT_SUCCESS);

int rval = EXIT_SUCCESS;
/*
* According Dicom standard (DICOM PS3.6 2016b - Data Dictionary)
* (0028, 0030) indicates physical X,Y spacing inside a slice;
* (0018, 0088) indicates physical Z spacing between slices;
* which above are also consistent with dcm2niix software.
* when we can not get (0018, 0088), we should compute spacing
* from the planes' positions (TODO, see PR 112).
* */
rval = this->GetElementDS<double>(0x0028, 0x0030, 2, _spacing, false);
if (rval != EXIT_SUCCESS)
{
Expand All @@ -1196,31 +1208,26 @@ DCMTKFileReader::GetSpacing(double * const spacing) const
rval = this->GetElementDS<double>(0x0018, 0x2010, 2, &_spacing[0], false);
}
}

if (rval == EXIT_SUCCESS)
{
// slice thickness
spacing[0] = _spacing[1];
spacing[1] = _spacing[0];
/*
* According Dicom standard (DICOM PS3.6 2016b - Data Dictionary)
* (0028, 0030) indicates physical X,Y spacing inside a slice;
* (0018, 0088) indicates physical Z spacing between slices;
* which above are also consistent with Dcom2iix software.
* when we can not get (0018, 0088), we should compute spacing
* from the planes' positions (TODO, see PR 112).
* */
if (GetElementDS<double>(0x0018, 0x0088, 1, &_spacing[2], false) == EXIT_SUCCESS)
{
spacing[2] = _spacing[2];
}
else
{
// punt, thicknes of 1
spacing[2] = 1.0;
}
return rval;
_have_val[0] = true;
_have_val[1] = true;
}

if (GetElementDS<double>(0x0018, 0x0088, 1, &_spacing[2], false) == EXIT_SUCCESS)
{
spacing[2] = _spacing[2];
_have_val[2] = true;
}

if (_have_val[0] && _have_val[1] && _have_val[2])
{
return EXIT_SUCCESS;
}

// this is for multiframe images -- preferentially use the shared
// functional group, and then the per-frame functional group
unsigned short candidateSequences[2] = {
Expand Down Expand Up @@ -1250,26 +1257,43 @@ DCMTKFileReader::GetSpacing(double * const spacing) const
* when we can not get (0018, 0088), we should compute spacing
* from the planes' positions (TODO, see PR 112).
* */
if (subSequence.GetElementDS<double>(0x0028, 0x0030, 2, _spacing, false) == EXIT_SUCCESS)
if (!(_have_val[0] && _have_val[1]))
{
spacing[0] = _spacing[1];
spacing[1] = _spacing[0];
if (subSequence.GetElementDS<double>(0x0018, 0x0088, 1, &_spacing[2], false) == EXIT_SUCCESS)
if (subSequence.GetElementDS<double>(0x0028, 0x0030, 2, _spacing, false) == EXIT_SUCCESS)
{
spacing[2] = _spacing[2];
spacing[0] = _spacing[1];
spacing[1] = _spacing[0];
_have_val[0] = true;
_have_val[1] = true;
}
else
}
if (!_have_val[2])
{
if (subSequence.GetElementDS<double>(0x0018, 0x0088, 1, &_spacing[2], false) == EXIT_SUCCESS)
{
// punt, zSpace of 1
spacing[2] = 1.0;
spacing[2] = _spacing[2];
_have_val[2] = true;
}
break;
}
if (_have_val[0] && _have_val[1] && _have_val[2])
{
return EXIT_SUCCESS;
}
}
}
}
}
return rval;

if (_have_val[0] && _have_val[1] && _have_val[2])
{
return EXIT_SUCCESS;
}

// not all (or no) spacing values could be found.
spacing[0] = _spacing[0];
spacing[1] = _spacing[1];
spacing[2] = _spacing[2];
return EXIT_FAILURE;
}

int
Expand Down
8 changes: 7 additions & 1 deletion Modules/IO/DCMTK/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,8 @@ set(ITKIODCMTKTests
itkDCMTKLoggerTest.cxx
itkDCMTKImageIOSlopeInterceptTest.cxx
itkDCMTKImageIOMultiFrameImageTest.cxx
itkDCMTKImageIONoPreambleTest.cxx)
itkDCMTKImageIONoPreambleTest.cxx
itkDCMTKImageIOSpacingTest.cxx)

createtestdriver(ITKIODCMTK "${ITKIODCMTK-Test_LIBRARIES}" "${ITKIODCMTKTests}")

Expand Down Expand Up @@ -188,3 +189,8 @@ itk_add_test(
ITKIODCMTKTestDriver
itkDCMTKImageIONoPreambleTest
DATA{Input/NoPreambleDicomTest.dcm})

itk_add_test(
NAME itkDCMTKImageIOSpacingTest
COMMAND ITKIODCMTKTestDriver itkDCMTKImageIOSpacingTest
DATA{Input/MultiFrameMRIZSpacing.dcm})
1 change: 1 addition & 0 deletions Modules/IO/DCMTK/test/Input/MultiFrameMRIZSpacing.dcm.cid
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
bafybeiabdmcpjhqoi7hfxx4vep27nc5y2ubclimgkcosk2cgpovbulwxim
172 changes: 172 additions & 0 deletions Modules/IO/DCMTK/test/itkDCMTKImageIOSpacingTest.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,172 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/

#include "itkImageFileReader.h"
#include "itkDCMTKImageIO.h"
#include "itkImageRegionConstIterator.h"
#include "itkMultiplyImageFilter.h"
#include "itkAddImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkStatisticsImageFilter.h"
#include "itkMath.h"
#include "itkTestingMacros.h"

using PixelType = short;
using ImageType = itk::Image<PixelType, 3>;
using DirectionType = ImageType::DirectionType;
using SpacingType = ImageType::SpacingType;
using PointType = ImageType::PointType;
using ReaderType = itk::ImageFileReader<ImageType>;
using ImageIOType = itk::DCMTKImageIO;

namespace
{
//
// if the difference is < average/10000, close enough
bool
CloseEnough(double a, double b)
{
double diff = itk::Math::abs(a - b);
double avg = (itk::Math::abs(a) + itk::Math::abs(b)) / 2.0;
if (diff == 0.0 || diff < avg / 10000.00)
{
return true;
}
return false;
}

bool
Equal(DirectionType dir1, DirectionType dir2)
{
for (unsigned int i = 0; i < 3; ++i)
{
for (unsigned int j = 0; j < 3; ++j)
{
if (!CloseEnough(dir1(i, j), dir2(i, j)))
{
return false;
}
}
}
return true;
}

bool
Equal(SpacingType spacing1, SpacingType spacing2)
{
for (unsigned int i = 0; i < 3; ++i)
{
if (!CloseEnough(spacing1[i], spacing2[i]))
{
return false;
}
}
return true;
}

bool
Equal(PointType p1, PointType p2)
{
for (unsigned int i = 0; i < 3; ++i)
{
if (!CloseEnough(p1[i], p2[i]))
{
return false;
}
}
return true;
}
} // namespace

// This tests if slice spacing is correctly read from a Philips MRI where
// SpacingBetweenSlices is stored as a global tag.
//
int
itkDCMTKImageIOSpacingTest(int argc, char * argv[])
{
if (argc < 2)
{
std::cerr << "Missing Parameters." << std::endl;
std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) << " multiframeImage" << std::endl;
return EXIT_FAILURE;
}

auto dcmImageIO = ImageIOType::New();

auto reader = ReaderType::New();
reader->SetFileName(argv[1]);
reader->SetImageIO(dcmImageIO);

ITK_TRY_EXPECT_NO_EXCEPTION(reader->Update());


const ImageType::Pointer im = reader->GetOutput();
const DirectionType dir = im->GetDirection();
const SpacingType spacing = im->GetSpacing();
const PointType origin = im->GetOrigin();

std::cerr << "Direction " << dir << std::endl
<< "Spacing " << spacing << std::endl
<< "Origin " << origin << std::endl;
DirectionType expectedDirection;
expectedDirection(0, 0) = 0.999894;
expectedDirection(0, 1) = 0.0145622;
expectedDirection(0, 2) = -0.000350456;
expectedDirection(1, 0) = 0.000350419;
expectedDirection(1, 1) = 5.1034e-06;
expectedDirection(1, 2) = 1.0;
expectedDirection(2, 0) = 0.0145622;
expectedDirection(2, 1) = -0.999894;
expectedDirection(2, 2) = -1.33323e-11;
if (!Equal(dir, expectedDirection))
{
std::cerr << "Expected directions" << std::endl
<< expectedDirection << std::endl
<< "Actual directions" << std::endl
<< dir << std::endl;
return EXIT_FAILURE;
}
SpacingType expectedSpacing;
expectedSpacing[0] = 0.55;
expectedSpacing[1] = 0.55;
expectedSpacing[2] = 1.3;
if (!Equal(spacing, expectedSpacing))
{
std::cerr << "Expected spacing" << std::endl
<< expectedSpacing << std::endl
<< "Actual spacing" << std::endl
<< spacing << std::endl;
return EXIT_FAILURE;
}
PointType expectedOrigin;
expectedOrigin[0] = -112.767;
expectedOrigin[1] = -41.4828;
expectedOrigin[2] = 117.567;
if (!Equal(origin, expectedOrigin))
{
std::cerr << "Expected origin" << std::endl
<< expectedOrigin << std::endl
<< "Actual origin" << std::endl
<< origin << std::endl;
return EXIT_FAILURE;
}


std::cout << "Test finished" << std::endl;
return EXIT_SUCCESS;
}

0 comments on commit dbd46b4

Please sign in to comment.