common/vvImageReader.cxx

   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 #ifndef VVIMAGEREADER_CXX
  19 #define VVIMAGEREADER_CXX
  20
  21 #include <itkImageFileReader.h>
  22 #include "vvImageReader.h"
  23 #include "vvImageReader.txx"
  24 #include "clitkTransformUtilities.h"
  25
  26 //------------------------------------------------------------------------------
  27 vvImageReader::vvImageReader()
  28 {
  29   mImage = NULL;
  30   mInputFilenames.resize(0);
  31   mLastError = "";
  32   mType = UNDEFINEDIMAGETYPE;
  33   mSlice = 0;
  34 }
  35 //------------------------------------------------------------------------------
  36
  37
  38 //------------------------------------------------------------------------------
  39 vvImageReader::~vvImageReader() { }
  40 //------------------------------------------------------------------------------
  41
  42
  43 //------------------------------------------------------------------------------
  44 void vvImageReader::Update()
  45 {
  46   Update(mType);
  47 }
  48 //------------------------------------------------------------------------------
  49
  50
  51 //------------------------------------------------------------------------------
  52 void vvImageReader::Update(LoadedImageType type)
  53 {
  54   itk::ImageIOBase::Pointer reader = itk::ImageIOFactory::CreateImageIO(mInputFilenames[0].c_str(), itk::ImageIOFactory::ReadMode);
  55   if (!reader) {
  56     mLastError="Unable to read file.";
  57   } else {
  58     reader->SetFileName(mInputFilenames[0]);
  59     reader->ReadImageInformation();
  60     if (mInputFilenames.size() > 1)
  61       Update(reader->GetNumberOfDimensions()+1,reader->GetComponentTypeAsString(reader->GetComponentType()),type);
  62     else
  63       Update(reader->GetNumberOfDimensions(),reader->GetComponentTypeAsString(reader->GetComponentType()),type);
  64   }
  65 }
  66 //------------------------------------------------------------------------------
  67
  68
  69 //------------------------------------------------------------------------------
  70 void vvImageReader::Update(int dim,std::string inputPixelType, LoadedImageType type)
  71 {
  72   //CALL_FOR_ALL_DIMS(dim,UpdateWithDim,inputPixelType);
  73   mType = type;
  74   mDim = dim;
  75   mInputPixelType=inputPixelType;
  76
  77   switch(mDim)     {
  78   case 2:
  79     UpdateWithDim<2>(mInputPixelType);
  80     break;
  81   case 3:
  82     UpdateWithDim<3>(mInputPixelType);
  83     break;
  84   case 4:
  85     UpdateWithDim<4>(mInputPixelType);
  86     break;
  87   default:
  88     std::cerr << "dimension unknown in Update ! " << std::endl;
  89   }
  90 }
  91 //------------------------------------------------------------------------------
  92
  93
  94 //------------------------------------------------------------------------------
  95 void vvImageReader::SetInputFilename(const std::string & filename)
  96 {
  97   mInputFilenames.resize(0);
  98   mInputFilenames.push_back(filename);
  99 }
 100 //------------------------------------------------------------------------------
 101
 102
 103 //------------------------------------------------------------------------------
 104 void vvImageReader::SetInputFilenames(const std::vector<std::string> & filenames)
 105 {
 106   mInputFilenames = filenames;
 107 }
 108 //------------------------------------------------------------------------------
 109
 110
 111 //------------------------------------------------------------------------------
 112 //Read transformation in NKI format (Xdr, transposed, cm)
 113 //void vvImageReader::ReadNkiImageTransform()
 114 //{
 115 //  bool bRead=false;
 116 //  std::string filename = mInputFilenames[0]+".MACHINEORIENTATION";
 117 //  if(itksys::SystemTools::FileExists(filename.c_str())){
 118 //    typedef itk::ImageFileReader< itk::Image< double, 2 > > MatrixReaderType;
 119 //    MatrixReaderType::Pointer readerTransfo = MatrixReaderType::New();
 120 //    readerTransfo->SetFileName(filename);
 121 //    try {
 122 //      bRead = true;
 123 //      readerTransfo->Update();
 124 //    } catch( itk::ExceptionObject & err ) {
 125 //      bRead=false;
 126 //      std::cerr << "Cannot read " << filename << std::endl
 127 //                << "The error is: " << err << std::endl;
 128 //    }
 129
 130 //    if (bRead) {
 131 //      //Transpose matrix (NKI format)
 132 //      for(int j=0; j<4; j++)
 133 //        for(int i=0; i<4; i++)
 134 //          mImage->GetTransform()->GetMatrix()->SetElement(j,i,readerTransfo->GetOutput()->GetBufferPointer()[4*i+j]);
 135
 136 //      //From cm to mm
 137 //      for(int i=0; i<3; i++)
 138 //        mImage->GetTransform()->GetMatrix()->SetElement(i,3,10*mImage->GetTransform()->GetMatrix()->GetElement(i,3));
 139
 140 //      mImage->GetTransform()->Inverse();
 141 //    }
 142 //  }
 143 //}
 144 //------------------------------------------------------------------------------
 145
 146
 147 //------------------------------------------------------------------------------
 148 //Read transformation in ASCII format
 149 void vvImageReader::ReadMatImageTransform()
 150 {
 151   std::string filename(mInputFilenames[0]);
 152   std::string ext(itksys::SystemTools::GetFilenameLastExtension(filename));
 153   if (ext.length() > 0) {
 154     size_t pos = filename.rfind(ext);
 155     filename.replace(pos, ext.length(), ".mat");
 156   }
 157   else
 158     filename += ".mat";
 159
 160   std::ifstream f(filename.c_str());
 161   if(f.is_open()) {
 162     f.close();
 163
 164     itk::Matrix<double, 4, 4> itkMat;
 165     itkMat.SetIdentity();
 166     try {
 167       itkMat = clitk::ReadMatrix3D(filename);
 168     }
 169     catch (itk::ExceptionObject & err) {
 170       itkWarningMacro(<< "Found " << filename
 171                       << " but this is not a 4x4 matrix so it is ignored.");
 172     }
 173
 174     vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
 175     matrix->Identity();
 176     for(int j=0; j<4; j++)
 177       for(int i=0; i<4; i++)
 178         matrix->SetElement(j,i,itkMat[j][i]);
 179
 180     // RP: 14/03/2011
 181     //  For some reason, the transformation matrix given in the MHD
 182     // file is inverted wrt the transformation given in the MAT file.
 183     // We don't really know where the inversion takes place inside
 184     // VTK or ITK. For what we could see in VV, the transformation
 185     // given in the MHD file seems "more correct" than that given in
 186     // the MAT file. For this reason, we invert the matrix read from
 187     // the MAT file before concatenating it to the current transformation.
 188     // Still, it would be nice to find out what happens exactly between
 189     // VTK and ITK...
 190     //
 191     // SR: 23/06/2011
 192     // We actually use vtkImageReslice which takes the inverse transform
 193     // as input. So it seems more correct to inverse the matrix given by
 194     // itk::ImageBase< VImageDimension >::GetDirection() than the matrix
 195     // in .mat which is indeed the inverse optimized by a typical
 196     // affine registration algorithm.
 197     // Validated using clitkAffineTransform --transform_grid
 198     if (matrix->Determinant() == 0)
 199     {
 200       vtkGenericWarningMacro("Matrix in " << filename.c_str() << " cannot be inverted (determinant = 0)");
 201     }
 202
 203     // TODO SR and BP: check on the list of transforms and not the first only
 204     mImage->GetTransform()[0]->PreMultiply();
 205     mImage->GetTransform()[0]->Concatenate(matrix);
 206     mImage->GetTransform()[0]->Update();
 207
 208     //for image sequences, apply the transform to each images of the sequence
 209     if (mImage->IsTimeSequence())
 210     {
 211         for (unsigned i = 1 ; i<mImage->GetTransform().size() ; i++)
 212         {
 213             mImage->GetTransform()[i]->PreMultiply();
 214             mImage->GetTransform()[i]->Concatenate(matrix);
 215             mImage->GetTransform()[i]->Update();
 216         }
 217     }
 218
 219   }
 220 }
 221 //------------------------------------------------------------------------------
 222 #endif
 223