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