]> Creatis software - clitk.git/blob - common/vvImageReader.cxx
Put back and refactored manual rigid registration
[clitk.git] / 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(itksys::SystemTools::GetFilenameWithoutExtension(mInputFilenames[0]));
152   filename += ".mat";
153
154   std::ifstream f(filename.c_str());
155   if(f.is_open()) {
156     f.close();
157
158     itk::Matrix<double, 4, 4> itkMat = clitk::ReadMatrix3D(filename);
159
160     vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
161     matrix->Identity();
162     for(int j=0; j<4; j++)
163       for(int i=0; i<4; i++)
164         matrix->SetElement(j,i,itkMat[j][i]);
165
166     // RP: 14/03/2011
167     //  For some reason, the transformation matrix given in the MHD
168     // file is inverted wrt the transformation given in the MAT file.
169     // We don't really know where the inversion takes place inside
170     // VTK or ITK. For what we could see in VV, the transformation
171     // given in the MHD file seems "more correct" than that given in
172     // the MAT file. For this reason, we invert the matrix read from
173     // the MAT file before concatenating it to the current transformation.
174     // Still, it would be nice to find out what happens exactly between
175     // VTK and ITK...
176     //
177     // SR: 23/06/2011
178     // We actually use vtkImageReslice which takes the inverse transform
179     // as input. So it seems more correct to inverse the matrix given by
180     // itk::ImageBase< VImageDimension >::GetDirection() than the matrix
181     // in .mat which is indeed the inverse optimized by a typical
182     // affine registration algorithm.
183     // Validated using clitkAffineTransform --transform_grid
184     if (matrix->Determinant() == 0)
185     {
186       vtkGenericWarningMacro("Matrix in " << filename.c_str() << " cannot be inverted (determinant = 0)");
187     }
188
189     mImage->GetTransform()->PostMultiply();
190     mImage->GetTransform()->Concatenate(matrix);
191     mImage->GetTransform()->Update();
192   }
193 }
194 //------------------------------------------------------------------------------
195 #endif
196