1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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
21 #include <itkImageFileReader.h>
22 #include "gdcmImageHelper.h"
23 #include "vvImageReader.h"
24 #include "vvImageReader.txx"
25 #include "clitkTransformUtilities.h"
26 #include "clitkElastix.h"
28 //------------------------------------------------------------------------------
29 vvImageReader::vvImageReader()
32 mInputFilenames.resize(0);
34 mType = UNDEFINEDIMAGETYPE;
36 mPatientCoordinateSystem = false;
38 //------------------------------------------------------------------------------
41 //------------------------------------------------------------------------------
42 vvImageReader::~vvImageReader() { }
43 //------------------------------------------------------------------------------
46 //------------------------------------------------------------------------------
47 void vvImageReader::Update()
51 //------------------------------------------------------------------------------
54 //------------------------------------------------------------------------------
55 void vvImageReader::SetPatientCoordinateSystem(bool patientCoordinateSystem)
57 mPatientCoordinateSystem = patientCoordinateSystem;
59 //------------------------------------------------------------------------------
62 //------------------------------------------------------------------------------
63 void vvImageReader::Update(LoadedImageType type)
65 itk::ImageIOBase::Pointer reader = itk::ImageIOFactory::CreateImageIO(mInputFilenames[0].c_str(), itk::ImageIOFactory::ReadMode);
67 mLastError="Unable to read file.";
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);
77 Update(reader->GetNumberOfDimensions(),reader->GetComponentTypeAsString(reader->GetComponentType()),type);
80 //------------------------------------------------------------------------------
83 //------------------------------------------------------------------------------
84 void vvImageReader::Update(int dim,std::string inputPixelType, LoadedImageType type)
86 //CALL_FOR_ALL_DIMS(dim,UpdateWithDim,inputPixelType);
89 mInputPixelType=inputPixelType;
93 UpdateWithDim<2>(mInputPixelType);
96 UpdateWithDim<3>(mInputPixelType);
99 UpdateWithDim<4>(mInputPixelType);
102 std::cerr << "dimension unknown in Update ! " << std::endl;
105 //------------------------------------------------------------------------------
108 //------------------------------------------------------------------------------
109 void vvImageReader::SetInputFilename(const std::string & filename)
111 mInputFilenames.resize(0);
112 mInputFilenames.push_back(filename);
114 //------------------------------------------------------------------------------
117 //------------------------------------------------------------------------------
118 void vvImageReader::SetInputFilenames(const std::vector<std::string> & filenames)
120 mInputFilenames = filenames;
122 //------------------------------------------------------------------------------
125 //------------------------------------------------------------------------------
126 //Read transformation in NKI format (Xdr, transposed, cm)
127 //void vvImageReader::ReadNkiImageTransform()
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);
137 // readerTransfo->Update();
138 // } catch( itk::ExceptionObject & err ) {
140 // std::cerr << "Cannot read " << filename << std::endl
141 // << "The error is: " << err << std::endl;
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]);
151 // for(int i=0; i<3; i++)
152 // mImage->GetTransform()->GetMatrix()->SetElement(i,3,10*mImage->GetTransform()->GetMatrix()->GetElement(i,3));
154 // mImage->GetTransform()->Inverse();
158 //------------------------------------------------------------------------------
161 //------------------------------------------------------------------------------
162 //Read transformation in ASCII format
163 void vvImageReader::ReadMatImageTransform()
165 std::string filename(mInputFilenames[0]);
166 std::string ext(itksys::SystemTools::GetFilenameLastExtension(filename));
168 // Try a ".mat" extension
169 if (ext.length() > 0) {
170 size_t pos = filename.rfind(ext);
171 filename.replace(pos, ext.length(), ".mat");
175 std::ifstream f(filename.c_str());
176 itk::Matrix<double, 4, 4> itkMat;
177 bool itkMatRead = false;
181 itkMat.SetIdentity();
183 itkMat = clitk::ReadMatrix3D(filename);
185 catch (itk::ExceptionObject & err) {
186 itkWarningMacro(<< "Found " << filename
187 << " but this is not a 4x4 matrix so it is ignored.");
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");
201 f.open(filename.c_str());
202 if(!itkMatRead && f.is_open()) {
204 itkMat = clitk::createMatrixFromElastixFile<3>(filename, true);
209 vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
211 for(int j=0; j<4; j++)
212 for(int i=0; i<4; i++)
213 matrix->SetElement(j,i,itkMat[j][i]);
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
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)
235 vtkGenericWarningMacro("Matrix in " << filename.c_str() << " cannot be inverted (determinant = 0)");
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();
243 //for image sequences, apply the transform to each images of the sequence
244 if (mImage->IsTimeSequence())
246 for (unsigned i = 1 ; i<mImage->GetTransform().size() ; i++)
248 mImage->GetTransform()[i]->PreMultiply();
249 mImage->GetTransform()[i]->Concatenate(matrix);
250 mImage->GetTransform()[i]->Update();
256 //------------------------------------------------------------------------------