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