]> Creatis software - clitk.git/blob - common/vvImageReader.txx
Add code to write dicom sequence tag in gateSimulation2Dicom
[clitk.git] / common / vvImageReader.txx
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
19 #ifndef VVIMAGEREADER_TXX
20 #define VVIMAGEREADER_TXX
21
22 #include <string>
23 #include <itkImageFileReader.h>
24 #include <itkImageSeriesReader.h>
25 #include <itkImageToVTKImageFilter.h>
26 #include <itkFlexibleVectorCastImageFilter.h>
27 #include "itkVectorImageToImageAdaptor.h"
28
29 #include <vtkTransform.h>
30
31 #include "clitkCommon.h"
32 #include "clitkConfiguration.h"
33 #include "vvFromITK.h"
34 //----------------------------------------------------------------------------
35 template<unsigned int VImageDimension>
36 void vvImageReader::UpdateWithDim(std::string InputPixelType)
37 {
38   if (mType == VECTORFIELD || mType == VECTORFIELDWITHTIME)
39   {
40     if (VImageDimension == 4)
41       UpdateWithDimAndInputVectorPixelType<itk::Vector<float,3>,VImageDimension>();
42     else
43       UpdateWithDimAndInputVectorPixelType<itk::Vector<float,VImageDimension>,VImageDimension>();
44   }
45   else if (InputPixelType == "short")
46     UpdateWithDimAndInputPixelType<short,VImageDimension>();
47   else if (InputPixelType == "unsigned_short")
48     UpdateWithDimAndInputPixelType<unsigned short,VImageDimension>();
49   else if (InputPixelType == "char")
50     UpdateWithDimAndInputPixelType<char,VImageDimension>();
51   else if (InputPixelType == "unsigned_char")
52     UpdateWithDimAndInputPixelType<unsigned char,VImageDimension>();
53   else if (InputPixelType == "int")
54     UpdateWithDimAndInputPixelType<int,VImageDimension>();
55   else if (InputPixelType == "unsigned_int")
56     UpdateWithDimAndInputPixelType<unsigned int,VImageDimension>();
57   else if (InputPixelType == "double")
58     UpdateWithDimAndInputPixelType<double,VImageDimension>();
59   else if (InputPixelType == "float")
60     UpdateWithDimAndInputPixelType<float,VImageDimension>();
61   else
62     std::cerr << "Error, input pixel type : " << InputPixelType << " unknown !" << std::endl;
63
64   if (mLastError.size()==0) {
65     //ReadNkiImageTransform();
66     ReadMatImageTransform();
67   }
68 }
69 //----------------------------------------------------------------------------
70
71
72 //----------------------------------------------------------------------------
73 template<class InputPixelType, unsigned int VImageDimension>
74 void vvImageReader::UpdateWithDimAndInputPixelType()
75 {
76   if (mType == MERGEDWITHTIME)   // In this case we can load the images
77     // one at the time to avoid excessive
78     // memory use
79   {
80     mImage=vvImage::New();
81
82     for (std::vector<std::string>::const_iterator i=mInputFilenames.begin(); i!=mInputFilenames.end(); i++) {
83       typedef itk::Image< InputPixelType, VImageDimension-1 > InputImageType;
84       typedef itk::ImageFileReader<InputImageType> ReaderType;
85       typename ReaderType::Pointer reader = ReaderType::New();
86       reader->ReleaseDataFlagOn();
87       reader->SetFileName(*i);
88       try {
89         mImage->AddItkImage<InputImageType>(reader->GetOutput());
90         mImage->ComputeScalarRangeBase<InputPixelType, VImageDimension-1>(reader->GetOutput());
91       } catch ( itk::ExceptionObject & err ) {
92         std::cerr << "Error while reading " << mInputFilenames[0].c_str()
93                   << " " << err << std::endl;
94         std::stringstream error;
95         error << err;
96         mLastError = error.str();
97         return;
98       }
99     }
100   } else if (mType == SLICED) {
101     mImage=vvImage::New();
102     typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
103     typedef itk::ImageFileReader<InputImageType> ReaderType;
104     typename ReaderType::Pointer reader = ReaderType::New();
105     reader->SetFileName(mInputFilenames[0]);
106     reader->UpdateOutputInformation();
107
108     typedef itk::Image< InputPixelType, VImageDimension-1 > SlicedImageType;
109     typedef itk::ExtractImageFilter<InputImageType, SlicedImageType> FilterType;
110
111     typename InputImageType::RegionType inputRegion = reader->GetOutput()->GetLargestPossibleRegion();
112     typename InputImageType::SizeType inputSize = inputRegion.GetSize();
113     typename InputImageType::IndexType start = inputRegion.GetIndex();
114     typename InputImageType::SizeType extractedRegionSize = inputSize;
115     typename InputImageType::RegionType extractedRegion;
116     extractedRegionSize[VImageDimension - 1] = 0;
117     extractedRegion.SetSize(extractedRegionSize);
118     start[VImageDimension - 1] = mSlice;
119     extractedRegion.SetIndex(start);
120
121     typename FilterType::Pointer filter = FilterType::New();
122     filter->SetExtractionRegion(extractedRegion);
123     filter->SetInput(reader->GetOutput());
124     filter->ReleaseDataFlagOn();
125     filter->SetDirectionCollapseToSubmatrix();
126     try {
127       mImage->AddItkImage<SlicedImageType>(filter->GetOutput());
128       mImage->ComputeScalarRangeBase<InputPixelType, VImageDimension-1>(filter->GetOutput());
129     } catch ( itk::ExceptionObject & err ) {
130       std::cerr << "Error while slicing " << mInputFilenames[0].c_str()
131                 << "(slice #" << mSlice << ") " << err << std::endl;
132       return;
133     }
134   } else if (mType == VECTORPIXELIMAGE) {
135     mImage=vvImage::New();
136     typedef itk::VectorImage< InputPixelType, VImageDimension-1 > InputImageType;
137     typedef itk::ImageFileReader<InputImageType> ReaderType;
138     typedef itk::Image<InputPixelType, VImageDimension> OutputImageType;
139     typename ReaderType::Pointer reader = ReaderType::New();
140     reader->SetFileName(mInputFilenames[0]);
141     reader->Update();
142     typename InputImageType::Pointer input= reader->GetOutput();
143
144     typedef itk::VectorImageToImageAdaptor<InputPixelType, VImageDimension-1> ImageAdaptorType;
145     typename ImageAdaptorType::Pointer adaptor = ImageAdaptorType::New();
146     typename OutputImageType::Pointer output = OutputImageType::New();
147
148     adaptor->SetExtractComponentIndex(0);
149     adaptor->SetImage(input);
150
151     //Create the output
152     typename OutputImageType::IndexType index;
153     index.Fill(0);
154     typename OutputImageType::SizeType size;
155     size.Fill(input->GetNumberOfComponentsPerPixel());
156     typename OutputImageType::SpacingType spacing;
157     spacing.Fill(1);
158     typename OutputImageType::PointType origin;
159     origin.Fill(0);
160     typename OutputImageType::DirectionType direction;
161     direction.SetIdentity();
162     for (unsigned int pixelDim=0; pixelDim<VImageDimension-1; ++pixelDim)
163     {
164       size[pixelDim]=adaptor->GetLargestPossibleRegion().GetSize(pixelDim);
165       spacing[pixelDim]=input->GetSpacing()[pixelDim];
166       origin[pixelDim]=input->GetOrigin()[pixelDim];
167       for (unsigned int pixelDim2=0; pixelDim2<VImageDimension-1; ++pixelDim2)
168       {
169         direction[pixelDim][pixelDim2]=input->GetDirection()[pixelDim][pixelDim2];
170       }
171     }
172     typename OutputImageType::RegionType region;
173     region.SetSize(size);
174     region.SetIndex(index);
175     output->SetRegions(region);
176     output->SetOrigin(origin);
177     output->SetSpacing(spacing);
178     output->SetDirection(direction);
179     output->Allocate();
180
181     //Copy each channel
182     for (unsigned int pixelDim=0; pixelDim<input->GetNumberOfComponentsPerPixel(); ++pixelDim)
183     {
184       adaptor->SetExtractComponentIndex(pixelDim);
185
186       itk::ImageRegionIterator<InputImageType> imageIterator(input,input->GetLargestPossibleRegion());
187
188       while(!imageIterator.IsAtEnd())
189       {
190         typename OutputImageType::IndexType indexVector;
191         indexVector.Fill(0);
192         for (unsigned int indexDim=0; indexDim<VImageDimension-1; ++indexDim)
193         {
194           indexVector[indexDim]=imageIterator.GetIndex().GetElement(indexDim);
195         }
196         indexVector[VImageDimension-1]=pixelDim;
197
198         output->SetPixel(indexVector, adaptor->GetPixel(imageIterator.GetIndex()));
199         ++imageIterator;
200       }
201     }
202
203 /*    if (VImageDimension == 4)
204       mType == VECTORPIXELIMAGEWITHTIME;
205     else
206       mType == VECTORPIXELIMAGE;*/
207
208     try {
209       mImage = vvImageFromITK<VImageDimension,InputPixelType>(output, mType == VECTORPIXELIMAGEWITHTIME);
210       mImage->ComputeScalarRangeBase<InputPixelType, VImageDimension>(output);
211     } catch ( itk::ExceptionObject & err ) {
212       std::cerr << "Error while slicing " << mInputFilenames[0].c_str()
213                 << " " << err << std::endl;
214       return;
215     }
216   } else {
217     if (mInputFilenames.size() > 1) {
218       typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
219       typedef itk::ImageSeriesReader<InputImageType> ReaderType;
220       typename ReaderType::Pointer reader = ReaderType::New();
221       reader->SetFileNames(mInputFilenames);
222       reader->ReleaseDataFlagOn();
223
224       try {
225         if (mType == IMAGEWITHTIME)
226         {
227           std::cerr << "We should never come here:" << std::endl
228             << "  Calling vvImageReader with multiple images and IMAGEWITHTIME is undefined." << std::endl
229             << "  You are probably looking for MERGEDWITHTIME Type." << std::endl;
230           return;
231         }
232         else
233           mImage=vvImageFromITK<VImageDimension,InputPixelType>(reader->GetOutput());
234       } catch ( itk::ExceptionObject & err ) {
235         std::cerr << "Error while reading image series:" << err << std::endl;
236         std::stringstream error;
237         error << err;
238         mLastError = error.str();
239         return;
240       }
241     } else {
242       typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
243       typedef itk::ImageFileReader<InputImageType> ReaderType;
244       typename ReaderType::Pointer reader = ReaderType::New();
245       reader->SetFileName(mInputFilenames[0]);
246       reader->ReleaseDataFlagOn();
247
248       try {
249         mImage = vvImageFromITK<VImageDimension,InputPixelType>(reader->GetOutput(), mType == IMAGEWITHTIME || mType == VECTORFIELDWITHTIME);
250       } catch ( itk::ExceptionObject & err ) {
251         std::cerr << "Error while reading " << mInputFilenames[0].c_str()
252                   << " " << err << std::endl;
253         std::stringstream error;
254         error << err;
255         mLastError = error.str();
256         return;
257       }
258     }
259   }
260 }
261 //----------------------------------------------------------------------------
262
263 //----------------------------------------------------------------------------
264 template<class InputPixelType, unsigned int VImageDimension>
265 void vvImageReader::UpdateWithDimAndInputVectorPixelType()
266 {
267   typedef itk::Image< InputPixelType, VImageDimension > InputImageType;
268   typename InputImageType::Pointer input;
269
270   if (mInputFilenames.size() > 1) {
271     typedef itk::ImageSeriesReader<InputImageType> ReaderType;
272     typename ReaderType::Pointer reader = ReaderType::New();
273     reader->SetFileNames(mInputFilenames);
274     reader->ReleaseDataFlagOn();
275     try {
276       reader->Update();
277       input = reader->GetOutput();
278     } catch ( itk::ExceptionObject & err ) {
279       std::cerr << "Error while reading image series:" << err << std::endl;
280       std::stringstream error;
281       error << err;
282       mLastError = error.str();
283       return;
284     }
285   } else {
286     typedef itk::ImageFileReader<InputImageType> ReaderType;
287     typename ReaderType::Pointer reader = ReaderType::New();
288     reader->SetFileName(mInputFilenames[0]);
289     reader->ReleaseDataFlagOn();
290     try {
291       reader->Update();
292       input = reader->GetOutput();
293     } catch ( itk::ExceptionObject & err ) {
294       std::cerr << "Error while reading " << mInputFilenames[0].c_str()
295         << " " << err << std::endl;
296       std::stringstream error;
297       error << err;
298       mLastError = error.str();
299       return;
300     }
301   }
302   
303   typedef itk::Image< itk::Vector<float , 3>, VImageDimension > VectorImageType;
304   typedef itk::FlexibleVectorCastImageFilter<InputImageType, VectorImageType> CasterType;
305   typename VectorImageType::Pointer casted_input;
306   typename CasterType::Pointer caster = CasterType::New();
307   caster->SetInput(input);
308   caster->Update();
309   casted_input = caster->GetOutput();
310   
311   mImage = vvImageFromITK<VImageDimension, itk::Vector<float, 3> >(casted_input, mType == IMAGEWITHTIME || mType == VECTORFIELDWITHTIME);
312 }
313 //----------------------------------------------------------------------------
314
315 #endif /* end #define vvImageReader_TXX */
316