X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=common%2FvvImageReader.cxx;h=2c094a8c3eda3cfb67e85acddafb9cceb1bf562c;hb=19b462067f2b839ebb7e7cc8ba84aa14f4152e8b;hp=0205355f6fafc40acd2dc82984ce11f1412dddc4;hpb=0fce6cfc8ec8c2177c876a9cac67048c9adc9a6c;p=clitk.git diff --git a/common/vvImageReader.cxx b/common/vvImageReader.cxx index 0205355..2c094a8 100644 --- a/common/vvImageReader.cxx +++ b/common/vvImageReader.cxx @@ -22,6 +22,7 @@ #include "vvImageReader.h" #include "vvImageReader.txx" #include "clitkTransformUtilities.h" +#include "clitkElastix.h" //------------------------------------------------------------------------------ vvImageReader::vvImageReader() @@ -77,13 +78,13 @@ void vvImageReader::Update(int dim,std::string inputPixelType, LoadedImageType t switch(mDim) { case 2: UpdateWithDim<2>(mInputPixelType); - break;; + break; case 3: UpdateWithDim<3>(mInputPixelType); - break;; + break; case 4: UpdateWithDim<4>(mInputPixelType); - break;; + break; default: std::cerr << "dimension unknown in Update ! " << std::endl; } @@ -110,37 +111,37 @@ void vvImageReader::SetInputFilenames(const std::vector & filenames //------------------------------------------------------------------------------ //Read transformation in NKI format (Xdr, transposed, cm) -void vvImageReader::ReadNkiImageTransform() -{ - bool bRead=false; - std::string filename = mInputFilenames[0]+".MACHINEORIENTATION"; - if(itksys::SystemTools::FileExists(filename.c_str())){ - typedef itk::ImageFileReader< itk::Image< double, 2 > > MatrixReaderType; - MatrixReaderType::Pointer readerTransfo = MatrixReaderType::New(); - readerTransfo->SetFileName(filename); - try { - bRead = true; - readerTransfo->Update(); - } catch( itk::ExceptionObject & err ) { - bRead=false; - std::cerr << "Cannot read " << filename << std::endl - << "The error is: " << err << std::endl; - } - - if (bRead) { - //Transpose matrix (NKI format) - for(int j=0; j<4; j++) - for(int i=0; i<4; i++) - mImage->GetTransform()->GetMatrix()->SetElement(j,i,readerTransfo->GetOutput()->GetBufferPointer()[4*i+j]); - - //From cm to mm - for(int i=0; i<3; i++) - mImage->GetTransform()->GetMatrix()->SetElement(i,3,10*mImage->GetTransform()->GetMatrix()->GetElement(i,3)); - - mImage->GetTransform()->Inverse(); - } - } -} +//void vvImageReader::ReadNkiImageTransform() +//{ +// bool bRead=false; +// std::string filename = mInputFilenames[0]+".MACHINEORIENTATION"; +// if(itksys::SystemTools::FileExists(filename.c_str())){ +// typedef itk::ImageFileReader< itk::Image< double, 2 > > MatrixReaderType; +// MatrixReaderType::Pointer readerTransfo = MatrixReaderType::New(); +// readerTransfo->SetFileName(filename); +// try { +// bRead = true; +// readerTransfo->Update(); +// } catch( itk::ExceptionObject & err ) { +// bRead=false; +// std::cerr << "Cannot read " << filename << std::endl +// << "The error is: " << err << std::endl; +// } + +// if (bRead) { +// //Transpose matrix (NKI format) +// for(int j=0; j<4; j++) +// for(int i=0; i<4; i++) +// mImage->GetTransform()->GetMatrix()->SetElement(j,i,readerTransfo->GetOutput()->GetBufferPointer()[4*i+j]); + +// //From cm to mm +// for(int i=0; i<3; i++) +// mImage->GetTransform()->GetMatrix()->SetElement(i,3,10*mImage->GetTransform()->GetMatrix()->GetElement(i,3)); + +// mImage->GetTransform()->Inverse(); +// } +// } +//} //------------------------------------------------------------------------------ @@ -150,19 +151,48 @@ void vvImageReader::ReadMatImageTransform() { std::string filename(mInputFilenames[0]); std::string ext(itksys::SystemTools::GetFilenameLastExtension(filename)); + + // Try a ".mat" extension if (ext.length() > 0) { size_t pos = filename.rfind(ext); filename.replace(pos, ext.length(), ".mat"); } else filename += ".mat"; - std::ifstream f(filename.c_str()); + itk::Matrix itkMat; + bool itkMatRead = false; if(f.is_open()) { - f.close(); + itkMatRead = true; + + itkMat.SetIdentity(); + try { + itkMat = clitk::ReadMatrix3D(filename); + } + catch (itk::ExceptionObject & err) { + itkWarningMacro(<< "Found " << filename + << " but this is not a 4x4 matrix so it is ignored."); + itkMatRead = false; + } + } + f.close(); - itk::Matrix itkMat = clitk::ReadMatrix3D(filename); + // Try a ".elx" extension + filename = mInputFilenames[0]; + if (ext.length() > 0) { + size_t pos = filename.rfind(ext); + filename.replace(pos, ext.length(), ".elx"); + } + else + filename += ".elx"; + f.open(filename.c_str()); + if(!itkMatRead && f.is_open()) { + itkMatRead = true; + itkMat = clitk::createMatrixFromElastixFile<3>(filename, true); + } + f.close(); + if(itkMatRead) { vtkSmartPointer matrix = vtkSmartPointer::New(); matrix->Identity(); for(int j=0; j<4; j++) @@ -192,9 +222,22 @@ void vvImageReader::ReadMatImageTransform() vtkGenericWarningMacro("Matrix in " << filename.c_str() << " cannot be inverted (determinant = 0)"); } - mImage->GetTransform()->PostMultiply(); - mImage->GetTransform()->Concatenate(matrix); - mImage->GetTransform()->Update(); + // TODO SR and BP: check on the list of transforms and not the first only + mImage->GetTransform()[0]->PreMultiply(); + mImage->GetTransform()[0]->Concatenate(matrix); + mImage->GetTransform()[0]->Update(); + + //for image sequences, apply the transform to each images of the sequence + if (mImage->IsTimeSequence()) + { + for (unsigned i = 1 ; iGetTransform().size() ; i++) + { + mImage->GetTransform()[i]->PreMultiply(); + mImage->GetTransform()[i]->Concatenate(matrix); + mImage->GetTransform()[i]->Update(); + } + } + } } //------------------------------------------------------------------------------