[Insight-users] Problem with DICOM writing

[Luis Ibanez]

Hi Markus,

Thanks for posting the additional details
and the screenshots from GLSliceViewer.

You are right, the problem seems to come
from the lack of updating the content of
the MetaDictionary.

Since you are resampling the image you must
update the values of:

      Origin
      Spacing
      Number of pixels


You will have an example on how to do this
in the code


    Insight/Examples/IO/
       DicomImageReadChangeHeaderWrite.cxx


This example is described in the ITK Software Guide


     http://www.itk.org/ItkSoftwareGuide.pdf



   Regards,


      Luis



---------------------------
Markus Weigert wrote:
>>
>> Hi Markus,
>>
>>
>> It doesn't seem that your problem is on the DICOM writing
>> but rather on the resampling.
>>
>>
>> Can you please save the resampled image using another
>> file format (not DICOM)  ?
>>
>> for example MetaImage, with the itk::ImageFileWriter...
>>
>> then look at the saved resampled image and check if it
>> looks the same as what you get with the DICOM writing.
>>
> 
> 
> Hi Luis,
> 
> thank you very much for your response.
> 
> When I write the image in Analyze - format for example,
> and read it again to render it with GLSliceView, everything is just fine.
> 
> I already figured that out.
> 
> 
> This is the reason, why I think, that the problem is in writing the 
> DICOM - serie.
> 
> I think the problem is, that I don't change the DICOM - tags  after 
> resampling the image
> before writing it out, instead I use the MetaDataDictionary associated 
> to the original
> FileSeriesReader (m_fileSeriesReader, see code below).
> 
> The resampling code was actually taken from one of your examples and 
> works very well
> (see attached screenshots), datatype used is short, internal datatype 
> used for resampling is float.
> 
> Hope this helps to figure out where the problem come's from.
> 
> 
> Thanks a lot,
> Markus
> 
> 
> 
> This is the code used for resampling:
> 
> typedef itk::Image<short,3> ImageType;
> 
> ImageType::Pointer RIPApplicationBase::resampleImage(ImageType::Pointer 
> inputImage,
>                double xFactor,
>                double yFactor,
>                double zFactor){
> 
> typedef itk::Image <float, 3>  InternalImageType;
> typedef itk::CastImageFilter<ImageType, InternalImageType> CastFilterType;
> typedef 
> itk::RecursiveGaussianImageFilter<InternalImageType,InternalImageType> 
> GaussianFilterType;
> 
> 
> CastFilterType::Pointer caster = CastFilterType::New();
> caster->SetInput(inputImage);
> 
> GaussianFilterType::Pointer smootherX = GaussianFilterType::New();
> GaussianFilterType::Pointer smootherY = GaussianFilterType::New();
> GaussianFilterType::Pointer smootherZ = GaussianFilterType::New();
> 
> smootherX->SetInput( caster->GetOutput() );
> smootherY->SetInput( smootherX->GetOutput() );
> smootherZ->SetInput( smootherY->GetOutput() );
> 
> 
> const ImageType::SpacingType& inputSpacing = inputImage->GetSpacing();
> const double sigmaX = inputSpacing[0] * xFactor;
> const double sigmaY = inputSpacing[1] * yFactor;
> const double sigmaZ = inputSpacing[2] * zFactor;
> 
> smootherX->SetSigma( sigmaX );
> smootherY->SetSigma( sigmaY );
> smootherZ->SetSigma( sigmaZ );
> smootherX->SetDirection( 0 );
> smootherY->SetDirection( 1 );
> smootherZ->SetDirection( 2 );
> smootherX->SetNormalizeAcrossScale( false );
> smootherY->SetNormalizeAcrossScale( false );
> smootherZ->SetNormalizeAcrossScale( false );
> 
> try
> {
>  smootherZ->Update();
> }
> catch( itk::ExceptionObject & excep )
> {
>  std::cerr << "Exception catched !" << std::endl;
>  std::cerr << excep << std::endl;
> }
> 
> std::cout << "Image Smoothed" << std::endl;
> 
> InternalImageType::ConstPointer smoothedImage = smootherZ->GetOutput();
> typedef itk::ResampleImageFilter<InternalImageType, ImageType > 
> ResampleFilterType;
> typedef itk::IdentityTransform<double, 3> TransformType;
> typedef itk::LinearInterpolateImageFunction<InternalImageType, double > 
> InterpolatorType;
> typedef ImageType::SizeType::SizeValueType SizeValueType;
> 
> 
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> resampler->SetInterpolator( interpolator );
> resampler->SetDefaultPixelValue(255); // value for regions without source
> ImageType::SpacingType spacing;
> spacing[0] = inputSpacing[0] * xFactor;
> spacing[1] = inputSpacing[1] * yFactor;
> spacing[2] = inputSpacing[2] * zFactor;
> resampler->SetOutputSpacing( spacing );
> 
> // Use the same origin
> 
> resampler->SetOutputOrigin( inputImage->GetOrigin() );
> ImageType::SizeType inputSize = 
> inputImage->GetLargestPossibleRegion().GetSize();
> 
> ImageType::SizeType size;
> size[0] = static_cast< SizeValueType >( inputSize[0] / xFactor );
> size[1] = static_cast< SizeValueType >( inputSize[1] / yFactor );
> size[2] = static_cast< SizeValueType >( inputSize[2] / zFactor );
> resampler->SetSize( size );
> resampler->SetInput( smoothedImage );
> 
> TransformType::Pointer transform = TransformType::New();
> transform->SetIdentity();
> resampler->SetTransform( transform );
> 
> resampler->Update();
> return resampler->GetOutput();
> 
> }
> 
> 
> 
> this is the code for writing the Dicom serie:
> 
> void RIPApplicationBase::writeMovingImage(std::string imageExtension, 
> std::string directory, std::string filename)
> {
> 
> typedef short PixelType;
> 
> //Write DICOM Series
> if(imageExtension =="dcm")
> {
> 
> typedef itk::Image<PixelType,2>       Image2DType;
> typedef itk::ImageSeriesWriter< ImageType, Image2DType > SeriesWriterType;
> typedef itk::NumericSeriesFileNames   NameGeneratorType;
> 
> std::string format  = filename + std::string("%03d") + std::string(".dcm");
> 
> SeriesWriterType::Pointer swriter = SeriesWriterType::New();
> movingImage->Update();
> swriter->SetInput( movingImage );
> 
> NameGeneratorType::Pointer nameGenerator = NameGeneratorType::New();
> ImageType::RegionType region = movingImage->GetLargestPossibleRegion();
> ImageType::IndexType start = region.GetIndex();
> ImageType::SizeType size = region.GetSize();
> const unsigned int firstSlice = start[2];
> const unsigned int lastSlice = start[2] + size[2] - 1;
> nameGenerator->SetStartIndex(firstSlice);
> nameGenerator->SetEndIndex(lastSlice);
> nameGenerator->SetIncrementIndex(1);
> 
>    nameGenerator->SetSeriesFormat(format.c_str());
>    swriter->SetFileNames( nameGenerator->GetFileNames()  );
> 
>  ImageIOType::Pointer gdcmIO = ImageIOType::New();
>  swriter->SetImageIO( gdcmIO );
>  swriter->SetMetaDataDictionaryArray( 
> m_fileSeriesReader->GetMetaDataDictionaryArray() );
> 
> // Try to write the serie:
>  try
>    {
>    swriter->Update();
>    }
>  catch( itk::ExceptionObject & excp )
>    {
>    std::cerr << "Exception thrown while writing the series " << std::endl;
>    std::cerr << excp << std::endl;
>    }
> 
> }
> 
> 
> 
> ----- Original Message ----- From: "Luis Ibanez" <luis.ibanez at kitware.com>
> To: "Markus Weigert" <m.weigert at fz-juelich.de>
> Cc: <insight-users at itk.org>
> Sent: Wednesday, September 28, 2005 8:05 PM
> Subject: Re: [Insight-users] Problem with DICOM writing
> 
> 
>>
>> Hi Markus,
>>
>>
>> It doesn't seem that your problem is on the DICOM writing
>> but rather on the resampling.
>>
>>
>> Can you please save the resampled image using another
>> file format (not DICOM)  ?
>>
>> for example MetaImage, with the itk::ImageFileWriter...
>>
>> then look at the saved resampled image and check if it
>> looks the same as what you get with the DICOM writing.
>>
>>
>> If so... then the problem is in your resampling code, and
>> we will need to look at your source code. Please post that
>> code to the list, and *PLEASE* make sure that you include
>> the declarations of the pixel types.
>>
>> [Posting a screenshot of what you see with GLSliceViewer
>>  will also help us to figure out what may be going wrong.]
>>
>>
>> A common reason of these two-images view is when there is
>> a mix of 8-bits and 16-bits pixels types.
>>
>>
>>
>>     Thanks
>>
>>
>>        Luis
>>
>>
>>
>> --------------------------
>> Markus Weigert wrote:
>>
>>> Hi insight-users,
>>>  I have a problem with writing Dicom files after reading
>>> a serie to an image and performing a resampling on the Image.
>>>  After the image is resampled, by decreasing resolution
>>> of a factor of 1/2 for example, I write the serie out.
>>>  However, when I read the written serie and render it with
>>> GLSliceView class e.g., it shows two small images side by side
>>> and the rest of the image is grey.
>>>  I may have missed to set some parameters for the output,
>>> but I can't determine which.
>>>  Can somebody please give me an advice, how to solve this problem.
>>>  Best regards,
>>> Markus
>>>  PS:
>>> Writing code is as follows:
>>>  //Write DICOM Series
>>>  if(imageExtension =="dcm")
>>>  {
>>>    typedef itk::Image<PixelType,2>       Image2DType;
>>>  typedef itk::ImageSeriesWriter< ImageType, Image2DType > 
>>> SeriesWriterType;
>>>  typedef itk::NumericSeriesFileNames   NameGeneratorType;
>>>  std::string format  = filename + std::string("%03d") + 
>>> std::string(".dcm");
>>>  SeriesWriterType::Pointer swriter = SeriesWriterType::New();
>>>  movingImage->Update();
>>>  swriter->SetInput( movingImage );
>>>  swriter->SetImageIO( m_gdcmIO );     //DICOM imageIO used for reading
>>>  NameGeneratorType::Pointer nameGenerator = NameGeneratorType::New();
>>>  ImageType::RegionType region = movingImage->GetLargestPossibleRegion();
>>>  ImageType::IndexType start = region.GetIndex();
>>>  ImageType::SizeType size = region.GetSize();
>>>  const unsigned int firstSlice = start[2];
>>>  const unsigned int lastSlice = start[2] + size[2] - 1;
>>>  nameGenerator->SetStartIndex(firstSlice);
>>>  nameGenerator->SetEndIndex(lastSlice);
>>>  nameGenerator->SetIncrementIndex(1);
>>>  nameGenerator->SetSeriesFormat(format.c_str());
>>>     swriter->SetFileNames( nameGenerator->GetFileNames()  );
>>>  m_seriesFileNames->SetOutputDirectory(directory.c_str());
>>>   swriter->SetMetaDataDictionaryArray( 
>>> m_fileSeriesReader->GetMetaDataDictionaryArray() ); 
>>> //FileSeriesReader which was used for reading
>>>   // Try to write the serie:
>>>   try
>>>     {
>>>     swriter->Update();
>>>     }
>>>   catch( itk::ExceptionObject & excp )
>>>     {
>>>     std::cerr << "Exception thrown while writing the series " << 
>>> std::endl;
>>>     std::cerr << excp << std::endl;
>>>     }
>>>  }
>>>
>>>
>>> ------------------------------------------------------------------------
>>>
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk.org
>>> http://www.itk.org/mailman/listinfo/insight-users
>>
>>
> 
> ------------------------------------------------------------------------
> 
> 
> ------------------------------------------------------------------------
> 
> 
> ------------------------------------------------------------------------
>