common/vvImageReader.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17   ===========================================================================**/
  18
  19 #ifndef VVIMAGEREADER_TXX
  20 #define VVIMAGEREADER_TXX
  21
  22 #include <string>
  23 #include <itkImageFileReader.h>
  24 #include <itkImageSeriesReader.h>
  25 #include <itkImageToVTKImageFilter.h>
  26 #include <itkFlexibleVectorCastImageFilter.h>
  27 #include "itkVectorImageToImageAdaptor.h"
  28
  29 #include <vtkTransform.h>
  30
  31 #include "clitkCommon.h"
  32 #include "clitkConfiguration.h"
  33 #include "vvFromITK.h"
  34 //----------------------------------------------------------------------------
  35 template<unsigned int VImageDimension>
  36 void vvImageReader::UpdateWithDim(std::string InputPixelType)
  37 {
  38   if (mType == VECTORFIELD || mType == VECTORFIELDWITHTIME)
  39   {
  40     if (VImageDimension == 4)
  41       UpdateWithDimAndInputVectorPixelType<itk::Vector<float,3>,VImageDimension>();
  42     else
  43       UpdateWithDimAndInputVectorPixelType<itk::Vector<float,VImageDimension>,VImageDimension>();
  44   }
  45   else if (InputPixelType == "short")
  46     UpdateWithDimAndInputPixelType<short,VImageDimension>();
  47   else if (InputPixelType == "unsigned_short")
  48     UpdateWithDimAndInputPixelType<unsigned short,VImageDimension>();
  49   else if (InputPixelType == "char")
  50     UpdateWithDimAndInputPixelType<char,VImageDimension>();
  51   else if (InputPixelType == "unsigned_char")
  52     UpdateWithDimAndInputPixelType<unsigned char,VImageDimension>();
  53   else if (InputPixelType == "int")
  54     UpdateWithDimAndInputPixelType<int,VImageDimension>();
  55   else if (InputPixelType == "unsigned_int")
  56     UpdateWithDimAndInputPixelType<unsigned int,VImageDimension>();
  57   else if (InputPixelType == "double")
  58     UpdateWithDimAndInputPixelType<double,VImageDimension>();
  59   else if (InputPixelType == "float")
  60     UpdateWithDimAndInputPixelType<float,VImageDimension>();
  61   else
  62     std::cerr << "Error, input pixel type : " << InputPixelType << " unknown !" << std::endl;
  63
  64   if (mLastError.size()==0) {
  65     //ReadNkiImageTransform();
  66     ReadMatImageTransform();
  67   }
  68 }
  69 //----------------------------------------------------------------------------
  70
  71
  72 //----------------------------------------------------------------------------
  73 template<class InputPixelType, unsigned int VImageDimension>
  74 void vvImageReader::UpdateWithDimAndInputPixelType()
  75 {
  76   if (mType == MERGEDWITHTIME)   // In this case we can load the images
  77     // one at the time to avoid excessive
  78     // memory use
  79   {
  80     mImage=vvImage::New();
  81
  82     for (std::vector<std::string>::const_iterator i=mInputFilenames.begin(); i!=mInputFilenames.end(); i++) {
  83       typedef itk::Image< InputPixelType, VImageDimension-1 > InputImageType;
  84       typedef itk::ImageFileReader<InputImageType> ReaderType;
  85       typename ReaderType::Pointer reader = ReaderType::New();
  86       reader->ReleaseDataFlagOn();
  87       reader->SetFileName(*i);
  88       try {
  89         mImage->AddItkImage<InputImageType>(reader->GetOutput());
  90         mImage->ComputeScalarRangeBase<InputPixelType, VImageDimension-1>(reader->GetOutput());
  91       } catch ( itk::ExceptionObject & err ) {
  92         std::cerr << "Error while reading " << mInputFilenames[0].c_str()
  93                   << " " << err << std::endl;
  94         std::stringstream error;
  95         error << err;
  96         mLastError = error.str();
  97         return;
  98       }
  99     }
 100   } else if (mType == SLICED) {
 101     mImage=vvImage::New();
 102     typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
 103     typedef itk::ImageFileReader<InputImageType> ReaderType;
 104     typename ReaderType::Pointer reader = ReaderType::New();
 105     reader->SetFileName(mInputFilenames[0]);
 106     reader->UpdateOutputInformation();
 107
 108     typedef itk::Image< InputPixelType, VImageDimension-1 > SlicedImageType;
 109     typedef itk::ExtractImageFilter<InputImageType, SlicedImageType> FilterType;
 110
 111     typename InputImageType::RegionType inputRegion = reader->GetOutput()->GetLargestPossibleRegion();
 112     typename InputImageType::SizeType inputSize = inputRegion.GetSize();
 113     typename InputImageType::IndexType start = inputRegion.GetIndex();
 114     typename InputImageType::SizeType extractedRegionSize = inputSize;
 115     typename InputImageType::RegionType extractedRegion;
 116     extractedRegionSize[VImageDimension - 1] = 0;
 117     extractedRegion.SetSize(extractedRegionSize);
 118     start[VImageDimension - 1] = mSlice;
 119     extractedRegion.SetIndex(start);
 120
 121     typename FilterType::Pointer filter = FilterType::New();
 122     filter->SetExtractionRegion(extractedRegion);
 123     filter->SetInput(reader->GetOutput());
 124     filter->ReleaseDataFlagOn();
 125     filter->SetDirectionCollapseToSubmatrix();
 126     try {
 127       mImage->AddItkImage<SlicedImageType>(filter->GetOutput());
 128       mImage->ComputeScalarRangeBase<InputPixelType, VImageDimension-1>(filter->GetOutput());
 129     } catch ( itk::ExceptionObject & err ) {
 130       std::cerr << "Error while slicing " << mInputFilenames[0].c_str()
 131                 << "(slice #" << mSlice << ") " << err << std::endl;
 132       return;
 133     }
 134   } else if (mType == VECTORPIXELIMAGE) {
 135     mImage=vvImage::New();
 136     typedef itk::VectorImage< InputPixelType, VImageDimension-1 > InputImageType;
 137     typedef itk::ImageFileReader<InputImageType> ReaderType;
 138     typedef itk::Image<InputPixelType, VImageDimension> OutputImageType;
 139     typename ReaderType::Pointer reader = ReaderType::New();
 140     reader->SetFileName(mInputFilenames[0]);
 141     reader->Update();
 142     typename InputImageType::Pointer input= reader->GetOutput();
 143
 144     typedef itk::VectorImageToImageAdaptor<InputPixelType, VImageDimension-1> ImageAdaptorType;
 145     typename ImageAdaptorType::Pointer adaptor = ImageAdaptorType::New();
 146     typename OutputImageType::Pointer output = OutputImageType::New();
 147
 148     adaptor->SetExtractComponentIndex(0);
 149     adaptor->SetImage(input);
 150
 151     //Create the output
 152     typename OutputImageType::IndexType index;
 153     index.Fill(0);
 154     typename OutputImageType::SizeType size;
 155     size.Fill(input->GetNumberOfComponentsPerPixel());
 156     typename OutputImageType::SpacingType spacing;
 157     spacing.Fill(1);
 158     typename OutputImageType::PointType origin;
 159     origin.Fill(0);
 160     typename OutputImageType::DirectionType direction;
 161     direction.SetIdentity();
 162     for (unsigned int pixelDim=0; pixelDim<VImageDimension-1; ++pixelDim)
 163     {
 164       size[pixelDim]=adaptor->GetLargestPossibleRegion().GetSize(pixelDim);
 165       spacing[pixelDim]=input->GetSpacing()[pixelDim];
 166       origin[pixelDim]=input->GetOrigin()[pixelDim];
 167       for (unsigned int pixelDim2=0; pixelDim2<VImageDimension-1; ++pixelDim2)
 168       {
 169         direction[pixelDim][pixelDim2]=input->GetDirection()[pixelDim][pixelDim2];
 170       }
 171     }
 172     typename OutputImageType::RegionType region;
 173     region.SetSize(size);
 174     region.SetIndex(index);
 175     output->SetRegions(region);
 176     output->SetOrigin(origin);
 177     output->SetSpacing(spacing);
 178     output->SetDirection(direction);
 179     output->Allocate();
 180
 181     //Copy each channel
 182     for (unsigned int pixelDim=0; pixelDim<input->GetNumberOfComponentsPerPixel(); ++pixelDim)
 183     {
 184       adaptor->SetExtractComponentIndex(pixelDim);
 185
 186       itk::ImageRegionIterator<InputImageType> imageIterator(input,input->GetLargestPossibleRegion());
 187
 188       while(!imageIterator.IsAtEnd())
 189       {
 190         typename OutputImageType::IndexType indexVector;
 191         indexVector.Fill(0);
 192         for (unsigned int indexDim=0; indexDim<VImageDimension-1; ++indexDim)
 193         {
 194           indexVector[indexDim]=imageIterator.GetIndex().GetElement(indexDim);
 195         }
 196         indexVector[VImageDimension-1]=pixelDim;
 197
 198         output->SetPixel(indexVector, adaptor->GetPixel(imageIterator.GetIndex()));
 199         ++imageIterator;
 200       }
 201     }
 202
 203 /*    if (VImageDimension == 4)
 204       mType == VECTORPIXELIMAGEWITHTIME;
 205     else
 206       mType == VECTORPIXELIMAGE;*/
 207
 208     try {
 209       mImage = vvImageFromITK<VImageDimension,InputPixelType>(output, mType == VECTORPIXELIMAGEWITHTIME);
 210       mImage->ComputeScalarRangeBase<InputPixelType, VImageDimension>(output);
 211     } catch ( itk::ExceptionObject & err ) {
 212       std::cerr << "Error while slicing " << mInputFilenames[0].c_str()
 213                 << " " << err << std::endl;
 214       return;
 215     }
 216   } else {
 217     if (mInputFilenames.size() > 1) {
 218       typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
 219       typedef itk::ImageSeriesReader<InputImageType> ReaderType;
 220       typename ReaderType::Pointer reader = ReaderType::New();
 221       reader->SetFileNames(mInputFilenames);
 222       reader->ReleaseDataFlagOn();
 223
 224       try {
 225         if (mType == IMAGEWITHTIME)
 226         {
 227           std::cerr << "We should never come here:" << std::endl
 228             << "  Calling vvImageReader with multiple images and IMAGEWITHTIME is undefined." << std::endl
 229             << "  You are probably looking for MERGEDWITHTIME Type." << std::endl;
 230           return;
 231         }
 232         else
 233           mImage=vvImageFromITK<VImageDimension,InputPixelType>(reader->GetOutput());
 234       } catch ( itk::ExceptionObject & err ) {
 235         std::cerr << "Error while reading image series:" << err << std::endl;
 236         std::stringstream error;
 237         error << err;
 238         mLastError = error.str();
 239         return;
 240       }
 241     } else {
 242       typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
 243       typedef itk::ImageFileReader<InputImageType> ReaderType;
 244       typename ReaderType::Pointer reader = ReaderType::New();
 245       reader->SetFileName(mInputFilenames[0]);
 246       reader->ReleaseDataFlagOn();
 247
 248       try {
 249         mImage = vvImageFromITK<VImageDimension,InputPixelType>(reader->GetOutput(), mType == IMAGEWITHTIME || mType == VECTORFIELDWITHTIME);
 250       } catch ( itk::ExceptionObject & err ) {
 251         std::cerr << "Error while reading " << mInputFilenames[0].c_str()
 252                   << " " << err << std::endl;
 253         std::stringstream error;
 254         error << err;
 255         mLastError = error.str();
 256         return;
 257       }
 258     }
 259   }
 260
 261   if (mType == DICOM && !mPatientCoordinateSystem)
 262     mImage->InitializeTransform();
 263 }
 264 //----------------------------------------------------------------------------
 265
 266 //----------------------------------------------------------------------------
 267 template<class InputPixelType, unsigned int VImageDimension>
 268 void vvImageReader::UpdateWithDimAndInputVectorPixelType()
 269 {
 270   typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
 271   typename InputImageType::Pointer input;
 272
 273   if (mInputFilenames.size() > 1) {
 274     typedef itk::ImageSeriesReader<InputImageType> ReaderType;
 275     typename ReaderType::Pointer reader = ReaderType::New();
 276     reader->SetFileNames(mInputFilenames);
 277     reader->ReleaseDataFlagOn();
 278     try {
 279       reader->Update();
 280       input = reader->GetOutput();
 281     } catch ( itk::ExceptionObject & err ) {
 282       std::cerr << "Error while reading image series:" << err << std::endl;
 283       std::stringstream error;
 284       error << err;
 285       mLastError = error.str();
 286       return;
 287     }
 288   } else {
 289     typedef itk::ImageFileReader<InputImageType> ReaderType;
 290     typename ReaderType::Pointer reader = ReaderType::New();
 291     reader->SetFileName(mInputFilenames[0]);
 292     reader->ReleaseDataFlagOn();
 293     try {
 294       reader->Update();
 295       input = reader->GetOutput();
 296     } catch ( itk::ExceptionObject & err ) {
 297       std::cerr << "Error while reading " << mInputFilenames[0].c_str()
 298         << " " << err << std::endl;
 299       std::stringstream error;
 300       error << err;
 301       mLastError = error.str();
 302       return;
 303     }
 304   }
 305
 306   typedef itk::Image< itk::Vector<float , 3>, VImageDimension > VectorImageType;
 307   typedef itk::FlexibleVectorCastImageFilter<InputImageType, VectorImageType> CasterType;
 308   typename VectorImageType::Pointer casted_input;
 309   typename CasterType::Pointer caster = CasterType::New();
 310   caster->SetInput(input);
 311   caster->Update();
 312   casted_input = caster->GetOutput();
 313
 314   mImage = vvImageFromITK<VImageDimension, itk::Vector<float, 3> >(casted_input, mType == IMAGEWITHTIME || mType == VECTORFIELDWITHTIME);
 315 }
 316 //----------------------------------------------------------------------------
 317
 318 #endif /* end #define vvImageReader_TXX */
 319