diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp index 0a01784044..2c9ef6df8b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp @@ -1,185 +1,222 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDiffusionDICOMFileReader.h" #include "mitkDiffusionDICOMFileReaderHelper.h" #include "mitkDiffusionHeaderSiemensDICOMFileReader.h" #include "mitkDiffusionHeaderSiemensMosaicDICOMFileReader.h" static void PerformHeaderAnalysis( mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType headers ) { unsigned int images = headers.size(); unsigned int unweighted_images = 0; unsigned int weighted_images = 0; mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType::const_iterator c_iter = headers.begin(); while( c_iter != headers.end() ) { const mitk::DiffusionImageDICOMHeaderInformation h = *c_iter; if( h.baseline ) unweighted_images++; if( h.b_value > 0 ) weighted_images++; ++c_iter; } MITK_INFO << " :: Analyzed volumes " << images << "\n" << " :: \t"<< unweighted_images << " b = 0" << "\n" << " :: \t"<< weighted_images << " b > 0"; } mitk::DiffusionDICOMFileReader::DiffusionDICOMFileReader() { } mitk::DiffusionDICOMFileReader::~DiffusionDICOMFileReader() { } bool mitk::DiffusionDICOMFileReader ::LoadImages() { // prepare data reading DiffusionDICOMFileReaderHelper helper; DiffusionDICOMFileReaderHelper::VolumeFileNamesContainer filenames; const size_t number_of_outputs = this->GetNumberOfOutputs(); for( size_t idx = 0; idx < number_of_outputs; idx++ ) { DICOMImageFrameList flist = this->GetOutput(idx).GetImageFrameList(); std::vector< std::string > FileNamesPerVolume; DICOMImageFrameList::const_iterator cIt = flist.begin(); while( cIt != flist.end() ) { FileNamesPerVolume.push_back( (*cIt)->Filename ); ++cIt; } filenames.push_back( FileNamesPerVolume ); } - helper.LoadToVector( filenames ); + // TODO : only prototyping to test loading of diffusion images + // we need some solution for the different types + typedef mitk::DiffusionImage DiffusionImageType; + DiffusionImageType::Pointer output_image = DiffusionImageType::New(); + DiffusionImageType::GradientDirectionContainerType::Pointer directions = + DiffusionImageType::GradientDirectionContainerType::New(); + double max_bvalue = 0; + for( size_t idx = 0; idx < number_of_outputs; idx++ ) + { + DiffusionImageDICOMHeaderInformation header = this->m_RetrievedHeader.at(idx); + + if( max_bvalue < header.b_value ) + max_bvalue = header.b_value; + } + + // normalize the retrieved gradient directions according to the set b-value (maximal one) + for( size_t idx = 0; idx < number_of_outputs; idx++ ) + { + DiffusionImageDICOMHeaderInformation header = this->m_RetrievedHeader.at(idx); + DiffusionImageType::GradientDirectionType grad = header.g_vector; + + grad.normalize(); + grad *= sqrt( header.b_value / max_bvalue ); + + directions->push_back( header.g_vector ); + } + + // initialize the output image + output_image->SetDirections( directions ); + output_image->SetB_Value( max_bvalue ); + output_image->SetVectorImage( helper.LoadToVector( filenames ) ); + output_image->InitializeFromVectorImage(); + output_image->UpdateBValueMap(); + + DICOMImageBlockDescriptor& block = this->InternalGetOutput(0); + + block.SetMitkImage( (mitk::Image::Pointer) output_image ); + } void mitk::DiffusionDICOMFileReader ::AnalyzeInputFiles() { Superclass::AnalyzeInputFiles(); // collect output from superclass size_t number_of_outputs = this->GetNumberOfOutputs(); if(number_of_outputs == 0) { MITK_ERROR << "Failed to parse input, retrieved 0 outputs from SeriesGDCMReader "; } DICOMImageBlockDescriptor block_0 = this->GetOutput(0); MITK_INFO << "Retrieved " << number_of_outputs << "outputs."; // collect vendor ID from the first output, first image StringList inputFilename; DICOMImageFrameInfo::Pointer frame_0 = block_0.GetImageFrameList().at(0); inputFilename.push_back( frame_0->Filename ); gdcm::Scanner gdcmScanner; gdcm::Tag t_vendor(0x008, 0x0070); gdcm::Tag t_imagetype(0x008, 0x008); // add DICOM Tag for vendor gdcmScanner.AddTag( t_vendor ); // add DICOM Tag for image type gdcmScanner.AddTag( t_imagetype ); gdcmScanner.Scan( inputFilename ); // retrieve both vendor and image type std::string vendor = gdcmScanner.GetValue( frame_0->Filename.c_str(), t_vendor ); std::string image_type = gdcmScanner.GetValue( frame_0->Filename.c_str(), t_imagetype ); MITK_INFO << "Got vendor: " << vendor << " image type " << image_type; mitk::DiffusionHeaderDICOMFileReader::Pointer header_reader; // parse vendor tag if( vendor.find("SIEMENS") != std::string::npos ) { if( image_type.find("MOSAIC") != std::string::npos ) { header_reader = mitk::DiffusionHeaderSiemensMosaicDICOMFileReader::New(); } else { header_reader = mitk::DiffusionHeaderSiemensDICOMFileReader::New(); } } else if( vendor.find("GE") != std::string::npos ) { } else if( vendor.find("PHILIPS") != std::string::npos ) { } else { // unknown vendor } if( header_reader.IsNull() ) { MITK_ERROR << "No header reader for given vendor. "; return; } bool canread = false; for( size_t idx; idx < number_of_outputs; idx++ ) { DICOMImageFrameInfo::Pointer frame = this->GetOutput( idx ).GetImageFrameList().at(0); canread = header_reader->ReadDiffusionHeader(frame->Filename); } // collect the information m_RetrievedHeader = header_reader->GetHeaderInformation(); // TODO : Analyze outputs + header information, i.e. for the loading confidence MITK_INFO << "----- Diffusion DICOM Analysis Report ---- "; PerformHeaderAnalysis( m_RetrievedHeader ); MITK_INFO << "==========================================="; } bool mitk::DiffusionDICOMFileReader ::CanHandleFile(const std::string &filename) { } diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h index b6861ee4af..be4c429bc1 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h @@ -1,101 +1,103 @@ #ifndef MITKDIFFUSIONDICOMFILEREADERHELPER_H #define MITKDIFFUSIONDICOMFILEREADERHELPER_H #include "mitkDiffusionImage.h" #include "itkImageSeriesReader.h" #include "itkVectorImage.h" namespace mitk { class DiffusionDICOMFileReaderHelper { public: typedef std::vector< std::string > StringContainer; typedef std::vector< StringContainer > VolumeFileNamesContainer; - - template< typename PixelType, unsigned int VecImageDimension> typename itk::VectorImage< PixelType, VecImageDimension >::Pointer LoadToVector( const VolumeFileNamesContainer& filenames //const itk::ImageBase<3>::RegionType requestedRegion ) { typedef itk::Image< PixelType, 3> InputImageType; typedef itk::ImageSeriesReader< InputImageType > SeriesReaderType; typename SeriesReaderType::Pointer probe_reader = SeriesReaderType::New(); probe_reader->SetFileNames( filenames.at(0) ); probe_reader->GenerateOutputInformation(); const itk::ImageBase<3>::RegionType requestedRegion = probe_reader->GetOutput()->GetLargestPossibleRegion(); MITK_INFO << " --- Probe reader : \n" << " Retrieved LPR " << requestedRegion; typedef itk::VectorImage< PixelType, 3 > VectorImageType; typename VectorImageType::Pointer output_image = VectorImageType::New(); output_image->SetNumberOfComponentsPerPixel( filenames.size() ); + output_image->SetSpacing( probe_reader->GetOutput()->GetSpacing() ); + output_image->SetOrigin( probe_reader->GetOutput()->GetOrigin() ); + output_image->SetDirection( probe_reader->GetOutput()->GetDirection() ); + output_image->SetLargestPossibleRegion( probe_reader->GetOutput()->GetLargestPossibleRegion() ); output_image->SetBufferedRegion( requestedRegion ); output_image->Allocate(); itk::ImageRegionIterator< VectorImageType > vecIter( output_image, requestedRegion ); VolumeFileNamesContainer::const_iterator volumesFileNamesIter = filenames.begin(); // iterate over the given volumes unsigned int component = 0; while( volumesFileNamesIter != filenames.end() ) { MITK_INFO << " ======== Loading volume " << component+1 << " of " << filenames.size(); typename SeriesReaderType::Pointer volume_reader = SeriesReaderType::New(); volume_reader->SetFileNames( *volumesFileNamesIter ); try { volume_reader->UpdateLargestPossibleRegion(); } catch( const itk::ExceptionObject &e) { mitkThrow() << " ITK Series reader failed : "<< e.what(); } itk::ImageRegionConstIterator< InputImageType > iRCIter ( volume_reader->GetOutput(), volume_reader->GetOutput()->GetLargestPossibleRegion() ); // transfer to vector image iRCIter.GoToBegin(); vecIter.GoToBegin(); while( !iRCIter.IsAtEnd() ) { typename VectorImageType::PixelType vector_pixel = vecIter.Get(); vector_pixel.SetElement( component, iRCIter.Get() ); vecIter.Set( vector_pixel ); ++vecIter; ++iRCIter; } ++volumesFileNamesIter; component++; } return output_image; } }; } #endif // MITKDIFFUSIONDICOMFILEREADERHELPER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp index 2a613f72b8..a6a37b1385 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp @@ -1,91 +1,110 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDiffusionDICOMFileReader.h" #include "mitkDiffusionDICOMFileReaderTestHelper.h" #include "mitkDICOMTagBasedSorter.h" #include "mitkDICOMSortByTag.h" +#include "mitkNrrdDiffusionImageWriter.h" + #include "mitkTestingMacros.h" using mitk::DICOMTag; int mitkDiffusionDICOMFileReaderTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkDiffusionDICOMFileReaderTest"); mitk::DiffusionDICOMFileReader::Pointer gdcmReader = mitk::DiffusionDICOMFileReader::New(); MITK_TEST_CONDITION_REQUIRED(gdcmReader.IsNotNull(), "DICOMITKSeriesGDCMReader can be instantiated."); mitk::DICOMFileReaderTestHelper::SetTestInputFilenames( argc,argv ); // check the Set/GetInput function mitk::DICOMFileReaderTestHelper::TestInputFilenames( gdcmReader ); MITK_INFO << "Test input filenanems"; // check that output is a good reproduction of input (no duplicates, no new elements) mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); MITK_INFO << "Test output"; // repeat test with some more realistic sorting gdcmReader = mitk::DiffusionDICOMFileReader::New(); // this also tests destruction mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New(); // Use tags as in Qmitk // all the things that split by tag in DicomSeriesReader tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0010) ); // Number of Rows tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0011) ); // Number of Columns tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0030) ); // Pixel Spacing tagSorter->AddDistinguishingTag( DICOMTag(0x0018, 0x1164) ); // Imager Pixel Spacing tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0037) ); // Image Orientation (Patient) // TODO add tolerance parameter (l. 1572 of original code) // TODO handle as real vectors! cluster with configurable errors! tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x000e) ); // Series Instance UID tagSorter->AddDistinguishingTag( DICOMTag(0x0018, 0x0050) ); // Slice Thickness tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0008) ); // Number of Frames tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0052) ); // Frame of Reference UID // gdcmReader->AddSortingElement( tagSorter ); //mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); mitk::DICOMSortCriterion::ConstPointer sorting = mitk::DICOMSortByTag::New( DICOMTag(0x0020, 0x0013), // instance number mitk::DICOMSortByTag::New( DICOMTag(0x0020, 0x0012), // aqcuisition number mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0032), // aqcuisition time mitk::DICOMSortByTag::New( DICOMTag(0x0018, 0x1060), // trigger time mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0018) // SOP instance UID (last resort, not really meaningful but decides clearly) ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer(); tagSorter->SetSortCriterion( sorting ); MITK_INFO << "Created sort"; gdcmReader->AddSortingElement( tagSorter ); mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); MITK_INFO << "Created sort"; //gdcmReader->PrintOutputs(std::cout, true); // really load images //mitk::DICOMFileReaderTestHelper::TestMitkImagesAreLoaded( gdcmReader ); gdcmReader->LoadImages(); + mitk::Image::Pointer loaded_image = gdcmReader->GetOutput(0).GetMitkImage(); + + mitk::DiffusionImage::Pointer d_img = static_cast*>( loaded_image.GetPointer() ); + + mitk::NrrdDiffusionImageWriter::Pointer writer = + mitk::NrrdDiffusionImageWriter::New(); + writer->SetFileName("/tmp/test.dwi"); + writer->SetInput(d_img ); + + try + { + writer->Update(); + } + catch( const itk::ExceptionObject& e) + { + MITK_TEST_FAILED_MSG( << "Writer failed : " << e.what() ); + } MITK_TEST_END(); }