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