[Insight-users] Image Intensities Changed During DICOM Read-Write-Read in ITK 3.20 vs. ITK 4.4.1
Bill Lorensen
bill.lorensen at gmail.com
Fri Aug 16 13:32:17 EDT 2013
I found the issue. It is more complicated than I first reported.
If you read, then write a dicom image (that has non trivial
slope/intercept) and write it without running any other ITK filters, the
writer writes the converted pixels but also writes the original
slope/intercept. Then, the reader reapplies the slope and intercept the
next time it is read.
So don't change the pixel type to float. I'll come up with a work around.
Eventually, I report and fix the bug.
Bill
On Fri, Aug 16, 2013 at 1:05 PM, Bill Lorensen <bill.lorensen at gmail.com>wrote:
> Your data has a non-zero intercept:
> (0028|1052) Rescale Intercept = 0
> (0028|1053) Rescale Slope = 3.010012
>
> When ITK reads DICOM images it applies the slope and intercept:
> value = slope * pixel + intercept
>
> So a short or unsigned short cannot represent the values exactly. But that
> should only have a small error if you use a signed quantity like short.
>
> If you define your image type as float, everything works correctly.
>
> ITK seems to read and convert the data, but is not writing it correctly.
>
> So, there does seem to be a problem with ITK, but at least for now you can
> declare your image type as float and all should be fine.
>
> Bill
>
>
> On Wed, Aug 14, 2013 at 5:22 AM, Constantine Zakkaroff <mnkz at leeds.ac.uk>wrote:
>
>> Hi Bill,
>>
>> I'm wondering if you've had a chance to try it out with the data I
>> emailed?
>>
>> I've tried it out with the 'Top to Bottom' dataset from
>> http://placid.nlm.nih.gov/**community/21<http://placid.nlm.nih.gov/community/21>and it all looks completely legit. But with my data it gets all screwed
>> up... What a strange thing!
>>
>> Regards,
>> Constantine
>>
>>
>>
>>
>> On 12/08/2013 20:09, Constantine Zakkaroff wrote:
>>
>>> Hi Bill.
>>>
>>> OK. I actually read and write only once slice in this test.
>>>
>>> Please find it attached: MR.1.3.46.670589.11.38102.5.0.**7452.rar.
>>>
>>> I get DICOM images without extension.
>>>
>>> I'm also attaching the saved file with changed intensity: file00.rar.
>>>
>>> Many thanks,
>>> Constantine
>>>
>>> On 12/08/2013 19:51, Bill Lorensen wrote:
>>>
>>>> I also tried VTK3.20 and got the same good results.
>>>>
>>>>
>>>> You could try one of these:
>>>> http://placid.nlm.nih.gov/**community/21<http://placid.nlm.nih.gov/community/21>
>>>>
>>>> If your data is anonymized, perhaps you can send me a link or just a few
>>>> slices.
>>>>
>>>> Bill
>>>>
>>>>
>>>>
>>>> On Mon, Aug 12, 2013 at 2:45 PM, Constantine Z <mnkz at leeds.ac.uk
>>>> <mailto:mnkz at leeds.ac.uk>> wrote:
>>>>
>>>> Thanks. Yes, I can see it all looks completely legit.
>>>>
>>>> Perhaps I can send you my data? MR from a Philips scanner. Or maybe
>>>> on Win it might work differently?
>>>>
>>>> Regards,
>>>> Constartine
>>>>
>>>>
>>>> On 12/08/2013 7:38 p.m., Bill Lorensen wrote:
>>>>
>>>>
>>>> I just tried your code on two different dicom datasets, a ct
>>>> and
>>>> an mr.
>>>>
>>>> Both worked for me.
>>>> ./DICOMTest ~/Data/itkTestcases/Back\ to\ Front/ .
>>>> ITK Version: 4.5
>>>> Stats for file: /home/lorensen/Data/__**itkTestcases/Back to
>>>> Front//i24389.MRDC.34
>>>> Min: 0, Max: 293, Mean: 28.7654, Sigma: 37.2513, Variance:
>>>> 1387.66, Sum:
>>>> 1.88517e+06
>>>> Stats for file: ./file00.dcm
>>>> Min: 0, Max: 293, Mean: 28.7654, Sigma: 37.2513, Variance:
>>>> 1387.66, Sum:
>>>> 1.88517e+06
>>>> [build] ./DICOMTest ~/Data/NAV/CT/ ./
>>>> ITK Version: 4.5
>>>> Stats for file: /home/lorensen/Data/NAV/CT//__**ct.1
>>>> Min: -1024, Max: 31744, Mean: 6277.24, Sigma: 13249.9,
>>>> Variance:
>>>> 1.75559e+08, Sum: 1.64554e+09
>>>> Stats for file: .//file00.dcm
>>>> Min: -1024, Max: 31744, Mean: 6277.24, Sigma: 13249.9,
>>>> Variance:
>>>> 1.75559e+08, Sum: 1.64554e+09
>>>>
>>>>
>>>>
>>>>
>>>> On Mon, Aug 12, 2013 at 11:46 AM, Constantine Zakkaroff
>>>> <mnkz at leeds.ac.uk <mailto:mnkz at leeds.ac.uk>
>>>> <mailto:mnkz at leeds.ac.uk <mailto:mnkz at leeds.ac.uk>>> wrote:
>>>>
>>>> Hi Brad,
>>>>
>>>> Yes, here's the minimal code:
>>>>
>>>>
>>>> ********************************____**************************
>>>> ****__**__************
>>>>
>>>> #include <cstdlib>
>>>> #include <itkImage.h>
>>>> #include <itkMetaImageIO.h>
>>>> #include <itkGDCMImageIO.h>
>>>> #include <itkGDCMSeriesFileNames.h>
>>>> #include <itkImageFileReader.h>
>>>> #include <itkImageFileWriter.h>
>>>> #include <itkStatisticsImageFilter.h>
>>>>
>>>> int main(int argc, char *argv[])
>>>> {
>>>> std::clog << "ITK Version: " <<
>>>> ITK_VERSION_STRING <<
>>>> std::endl;
>>>>
>>>> typedef itk::GDCMSeriesFileNames
>>>> NamesGeneratorType;
>>>> typedef itk::Image<short, 3> ImageType;
>>>>
>>>> NamesGeneratorType::Pointer namesGenerator =
>>>> NamesGeneratorType::New();
>>>> namesGenerator->____**SetUseSeriesDetails(true);
>>>> namesGenerator->____**
>>>> SetInputDirectory(argv[1]);
>>>>
>>>> typedef itk::StatisticsImageFilter<___**
>>>> _ImageType>
>>>> StatsFilterType;
>>>> typedef itk::ImageFileReader<____**ImageType>
>>>> ImageReaderType;
>>>> typedef itk::ImageFileWriter<____**ImageType>
>>>> ImageWriterType;
>>>>
>>>>
>>>> ImageReaderType::Pointer imageReader0 =
>>>> ImageReaderType::New();
>>>> imageReader0->SetImageIO(itk::**
>>>> ____GDCMImageIO::New());
>>>>
>>>> StatsFilterType::Pointer statsFilter0 =
>>>> StatsFilterType::New();
>>>>
>>>>
>>>> imageReader0->SetFileName(____**namesGenerator->____**
>>>> GetInputFileNames()[0]);
>>>> statsFilter0->SetInput(____**
>>>> imageReader0->GetOutput());
>>>>
>>>> statsFilter0->Update();
>>>>
>>>> std::clog << "Stats for file: " <<
>>>> namesGenerator->____**GetInputFileNames()[0] <<
>>>> std::endl;
>>>>
>>>> std::clog <<
>>>> "Min: " << statsFilter0->GetMinimum() <<
>>>> ", "
>>>> "Max: " << statsFilter0->GetMaximum() <<
>>>> ", "
>>>> "Mean: " << statsFilter0->GetMean() <<
>>>> ", "
>>>> "Sigma: " << statsFilter0->GetSigma() <<
>>>> ", "
>>>> "Variance: " <<
>>>> statsFilter0->GetVariance() << ", "
>>>> "Sum: " << statsFilter0->GetSum()
>>>> << std::endl;
>>>>
>>>> ImageWriterType::Pointer imageWriter0 =
>>>> ImageWriterType::New();
>>>> imageWriter0->SetInput(____**
>>>> imageReader0->GetOutput());
>>>> imageWriter0->SetImageIO(itk::**
>>>> ____GDCMImageIO::New());
>>>>
>>>> std::string fileName0 = std::string(argv[2]) +
>>>> "/file00.dcm";
>>>> imageWriter0->SetFileName(____**fileName0);
>>>>
>>>> imageWriter0->Update();
>>>>
>>>> ImageReaderType::Pointer imageReader1 =
>>>> ImageReaderType::New();
>>>> imageReader1->SetImageIO(itk::**
>>>> ____GDCMImageIO::New());
>>>>
>>>> StatsFilterType::Pointer statsFilter1 =
>>>> StatsFilterType::New();
>>>> imageReader1->SetFileName(____**fileName0);
>>>> statsFilter1->SetInput(____**
>>>> imageReader1->GetOutput());
>>>>
>>>> statsFilter1->Update();
>>>>
>>>> std::clog << "Stats for file: " << fileName0 <<
>>>> std::endl;
>>>> std::clog <<
>>>> "Min: " << statsFilter1->GetMinimum() <<
>>>> ", "
>>>> "Max: " << statsFilter1->GetMaximum() <<
>>>> ", "
>>>> "Mean: " << statsFilter1->GetMean() <<
>>>> ", "
>>>> "Sigma: " << statsFilter1->GetSigma() <<
>>>> ", "
>>>> "Variance: " <<
>>>> statsFilter1->GetVariance() << ", "
>>>> "Sum: " << statsFilter1->GetSum()
>>>> << std::endl;
>>>>
>>>> system("pause");
>>>>
>>>> return EXIT_SUCCESS;
>>>> }
>>>>
>>>> ********************************____**************************
>>>> ****__**__************
>>>>
>>>>
>>>> Many thanks,
>>>> Constantine
>>>>
>>>>
>>>> On 12/08/2013 14:48, Bradley Lowekamp wrote:
>>>>
>>>>
>>>> Do you have the code to share so that it can be easily
>>>> reproduced?
>>>>
>>>>
>>>> Thanks,
>>>> Brad
>>>>
>>>> On Aug 12, 2013, at 9:33 AM, Constantine Zakkaroff
>>>> <mnkz at leeds.ac.uk <mailto:mnkz at leeds.ac.uk>
>>>> <mailto:mnkz at leeds.ac.uk <mailto:mnkz at leeds.ac.uk>>> wrote:
>>>>
>>>> Hello ALL,
>>>>
>>>> I'm wondering if anyone can explain what's going
>>>> on
>>>> with the
>>>> changing image intensities in a read-write-read
>>>> round-trip.
>>>> I think it's a bug, but the issues.itk.org
>>>> <http://issues.itk.org>
>>>> <http://issues.itk.org> is unresponsive at the
>>>> moment.
>>>>
>>>> Another bug? ;)
>>>>
>>>> Anyhow, the strange thing is that the results for
>>>> the round
>>>> trip read-write-read test are very different for
>>>> ITK 3.20
>>>> and ITK 4.4.1.
>>>>
>>>> I use to StatisticsImageFilter to examine the
>>>> intensity
>>>> values. Here's what happens, when I simply read a
>>>> DICOM
>>>> image, print stats, save the image and read it
>>>> again,
>>>> followed by printing the stats:
>>>>
>>>> ITK Version: 3.20
>>>> Stats for original file:
>>>> Min: 0, Max: 2474, Mean: 169.236, Sigma: 256.719,
>>>> Variance:
>>>> 65904.7, Sum: 1.40371e+007
>>>> Stats for saved file:
>>>> Min: 0, Max: 2471, Mean: 167.132, Sigma: 256.125,
>>>> Variance:
>>>> 65600.1, Sum: 1.38626e+007
>>>>
>>>> ITK Version: 4.4
>>>> Stats for original file:
>>>> Min: 0, Max: 2474, Mean: 169.236, Sigma: 256.719,
>>>> Variance:
>>>> 65904.7, Sum: 1.40371e+007
>>>> Stats for saved file:
>>>> Min: 0, Max: 7446, Mean: 509.114, Sigma: 772.59,
>>>> Variance:
>>>> 596895, Sum: 4.22279e+007
>>>>
>>>> In both cases the code is identical. Just compiled
>>>> against
>>>> different versions of ITK. The images are read and
>>>> written
>>>> as itk::Image<unsigned short, 3>, with no
>>>> processing, no
>>>> changes to any aspect of the images.
>>>>
>>>> So, ITK 3.20 introduces a small change, which I'm
>>>> pretty
>>>> sure shouldn't be there. But in ITK 4.4.1 the
>>>> scale of
>>>> intensities changes three-fold.
>>>>
>>>> I'd not be surprised by a change of this magnitude
>>>> if I was
>>>> photocopying those images, but not when using
>>>> digital
>>>> software and hardware. :^) Is there something I'm
>>>> not
>>>> understanding here?
>>>>
>>>> None of this happens with
>>>> read(DICOM)-write(MHD)-read(__**__MHD), where the
>>>> original and
>>>>
>>>> round-trip values are identical:
>>>> Stats for saved file: file0.mhd
>>>> Min: 0, Max: 2474, Mean: 169.236, Sigma: 256.719,
>>>> Variance:
>>>> 65904.7, Sum: 1.40371e+007
>>>>
>>>> Many thanks,
>>>> Constantine
>>>>
>>>> ______________________________**___________
>>>>
>>>> Powered by www.kitware.com <
>>>> http://www.kitware.com>
>>>> <http://www.kitware.com>
>>>>
>>>> Visit other Kitware open-source projects at
>>>> http://www.kitware.com/____**opensource/opensource.html<http://www.kitware.com/____opensource/opensource.html>
>>>> <http://www.kitware.com/__**opensource/opensource.html<http://www.kitware.com/__opensource/opensource.html>
>>>> >
>>>>
>>>>
>>>> <http://www.kitware.com/__**opensource/opensource.html<http://www.kitware.com/__opensource/opensource.html>
>>>> <http://www.kitware.com/**opensource/opensource.html<http://www.kitware.com/opensource/opensource.html>
>>>> >>
>>>>
>>>> Kitware offers ITK Training Courses, for more
>>>> information visit:
>>>> http://www.kitware.com/____**products/protraining.php<http://www.kitware.com/____products/protraining.php>
>>>> <http://www.kitware.com/__**products/protraining.php<http://www.kitware.com/__products/protraining.php>
>>>> >
>>>>
>>>> <http://www.kitware.com/__**
>>>> products/protraining.php<http://www.kitware.com/__products/protraining.php>
>>>> <http://www.kitware.com/**products/protraining.php<http://www.kitware.com/products/protraining.php>
>>>> >>
>>>>
>>>> Please keep messages on-topic and check the ITK
>>>> FAQ at:
>>>> http://www.itk.org/Wiki/ITK___**__FAQ<http://www.itk.org/Wiki/ITK_____FAQ>
>>>> <http://www.itk.org/Wiki/ITK__**_FAQ<http://www.itk.org/Wiki/ITK___FAQ>
>>>> >
>>>>
>>>> <http://www.itk.org/Wiki/ITK__**_FAQ<http://www.itk.org/Wiki/ITK___FAQ>
>>>> <http://www.itk.org/Wiki/ITK_**FAQ<http://www.itk.org/Wiki/ITK_FAQ>
>>>> >>
>>>>
>>>> Follow this link to subscribe/unsubscribe:
>>>> http://www.itk.org/mailman/___**_listinfo/insight-users<http://www.itk.org/mailman/____listinfo/insight-users>
>>>> <http://www.itk.org/mailman/__**listinfo/insight-users<http://www.itk.org/mailman/__listinfo/insight-users>
>>>> >
>>>>
>>>> <http://www.itk.org/mailman/__**listinfo/insight-users<http://www.itk.org/mailman/__listinfo/insight-users>
>>>> <http://www.itk.org/mailman/**listinfo/insight-users<http://www.itk.org/mailman/listinfo/insight-users>
>>>> >>
>>>>
>>>>
>>>> ______________________________**___________
>>>>
>>>> Powered by www.kitware.com <http://www.kitware.com>
>>>> <http://www.kitware.com>
>>>>
>>>> Visit other Kitware open-source projects at
>>>> http://www.kitware.com/____**opensource/opensource.html<http://www.kitware.com/____opensource/opensource.html>
>>>> <http://www.kitware.com/__**opensource/opensource.html<http://www.kitware.com/__opensource/opensource.html>
>>>> >
>>>>
>>>> <http://www.kitware.com/__**opensource/opensource.html<http://www.kitware.com/__opensource/opensource.html>
>>>> <http://www.kitware.com/**opensource/opensource.html<http://www.kitware.com/opensource/opensource.html>
>>>> >>
>>>>
>>>> Kitware offers ITK Training Courses, for more information
>>>> visit:
>>>> http://www.kitware.com/____**products/protraining.php<http://www.kitware.com/____products/protraining.php>
>>>> <http://www.kitware.com/__**products/protraining.php<http://www.kitware.com/__products/protraining.php>
>>>> >
>>>>
>>>> <http://www.kitware.com/__**products/protraining.php<http://www.kitware.com/__products/protraining.php>
>>>> <http://www.kitware.com/**products/protraining.php<http://www.kitware.com/products/protraining.php>
>>>> >>
>>>>
>>>> Please keep messages on-topic and check the ITK FAQ at:
>>>> http://www.itk.org/Wiki/ITK___**__FAQ<http://www.itk.org/Wiki/ITK_____FAQ>
>>>> <http://www.itk.org/Wiki/ITK__**_FAQ<http://www.itk.org/Wiki/ITK___FAQ>
>>>> >
>>>> <http://www.itk.org/Wiki/ITK__**_FAQ<http://www.itk.org/Wiki/ITK___FAQ>
>>>> <http://www.itk.org/Wiki/ITK_**FAQ<http://www.itk.org/Wiki/ITK_FAQ>
>>>> >>
>>>>
>>>>
>>>> Follow this link to subscribe/unsubscribe:
>>>> http://www.itk.org/mailman/___**_listinfo/insight-users<http://www.itk.org/mailman/____listinfo/insight-users>
>>>> <http://www.itk.org/mailman/__**listinfo/insight-users<http://www.itk.org/mailman/__listinfo/insight-users>
>>>> >
>>>>
>>>> <http://www.itk.org/mailman/__**listinfo/insight-users<http://www.itk.org/mailman/__listinfo/insight-users>
>>>> <http://www.itk.org/mailman/**listinfo/insight-users<http://www.itk.org/mailman/listinfo/insight-users>
>>>> >>
>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> Unpaid intern in BillsBasement at noware dot com
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> Unpaid intern in BillsBasement at noware dot com
>>>>
>>>
>
>
> --
> Unpaid intern in BillsBasement at noware dot com
>
--
Unpaid intern in BillsBasement at noware dot com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20130816/38b6dc15/attachment.htm>
More information about the Insight-users
mailing list