]> Creatis software - clitk.git/blob - tools/clitkDicom2Image.cxx
Toggle interpolation for fusion and overlay actors to (when applicable)
[clitk.git] / tools / clitkDicom2Image.cxx
1 /*=========================================================================
2                                                                                 
3 Program:   clitk
4 Module:    $RCSfile: clitkDicom2Image.cxx,v $
5 Language:  C++
6 Date:      $Date: 2010/11/05 00:05:00 $
7 Version:   $Revision: 1.1 $
8                                                                                 
9 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
10 l'Image). All rights reserved. See Doc/License.txt or
11 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
12                                                                                 
13 This software is distributed WITHOUT ANY WARRANTY; without even
14 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 PURPOSE.  See the above copyright notices for more information.
16                                                                              
17 =========================================================================*/
18
19 /**
20    =================================================
21    * @file   clitkDicom2Image.cxx
22    * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
23    * @date   26 Oct 2006
24    * 
25    * @brief  
26    * 
27    =================================================*/
28
29 // clitk includes
30 #include "clitkDicom2Image_ggo.h"
31 #include <clitkCommon.h>
32 #include "clitkImageCommon.h"
33
34 // itk include
35 #include "itkImageRegionIterator.h"
36 #include "itkImageRegionConstIterator.h"
37
38 // Why this is needed ???
39 //const double itk::NumericTraits<double>::Zero = 0.0;
40
41 //====================================================================
42 int main(int argc, char * argv[]) {
43
44   // init command line
45   GGO(clitkDicom2Image, args_info);
46   std::vector<std::string> input_files;
47   ///if std_input is given, read the input files from stdin
48   if (args_info.std_input_given)
49   {
50       while (true)
51       {
52           std::string tmp;
53           std::cin >> tmp;
54           if(std::cin.good())
55               input_files.push_back(tmp);
56           else break;
57       }
58       args_info.inputs_num=input_files.size();
59   }
60   else for (unsigned int i=0;i<args_info.inputs_num;i++)
61       input_files.push_back(args_info.inputs[i]);
62
63   
64   //===========================================
65   /// Get slices locations ...
66     std::vector<double> sliceLocations;
67     for(unsigned int i=0; i<args_info.inputs_num; i++) {
68       //std::cout << "Reading <" << input_files[i] << std::endl;        
69       gdcm::File * header = clitk::readDicomHeader(input_files[i]);
70       sliceLocations.push_back(header->GetZOrigin());
71       // check
72       /*if (!header->IsSignedPixelData()) {
73         std::cerr << "Pixel type not signed ! " << std::endl;
74         std::cerr << "In file " << input_files[i] << std::endl;
75         exit(0);
76         }*/
77       if (header->GetPixelSize() != 2) {
78         std::cerr << "Pixel type 2 bytes ! " << std::endl;
79         std::cerr << "In file " << input_files[i] << std::endl;
80         exit(0);
81       }
82     }
83
84     //===========================================
85     // Sort slices locations ...
86     std::vector<int> sliceIndex;
87     clitk::GetSortedIndex(sliceLocations, sliceIndex);
88     if (args_info.verboseSliceLocation_flag) {
89       std::cout << sliceLocations[sliceIndex[0]] << " -> " 
90                 << sliceIndex[0] << " / " << 0 << " => " 
91                 << "0 mm " 
92                 << input_files[sliceIndex[0]]
93                 << std::endl;
94       for(unsigned int i=1; i<sliceIndex.size(); i++) {
95         std::cout << sliceLocations[sliceIndex[i]] << " -> " 
96                   << sliceIndex[i] << " / " << i << " => " 
97                   << sliceLocations[sliceIndex[i]] - sliceLocations[sliceIndex[i-1]] 
98                   << "mm " 
99                   << input_files[sliceIndex[i]]
100                   << std::endl;
101       }
102     }
103
104     //===========================================
105     // Analyze slices locations ...
106     double currentDist;
107     double dist=0;
108     double tolerance = args_info.tolerance_arg;
109     double previous = sliceLocations[sliceIndex[0]];
110     for(unsigned int i=1; i<sliceIndex.size(); i++) {
111       currentDist = sliceLocations[sliceIndex[i]]-previous;
112       if (i!=1) {
113         if (fabs(dist-currentDist) > tolerance) {
114           std::cout << "ERROR : " << std::endl
115                     << "Current slice pos is  = " << sliceLocations[sliceIndex[i]] << std::endl
116                     << "Previous slice pos is = " << previous << std::endl
117                     << "Current file is       = " << input_files[sliceIndex[i]] << std::endl
118                     << "Current index is      = " << i << std::endl
119                     << "Current sortindex is  = " << sliceIndex[i] << std::endl
120                     << "Current slice diff    = " << dist << std::endl
121                     << "Current error         = " << fabs(dist-currentDist) << std::endl;
122           exit(1);
123         }
124       }
125       else dist = currentDist;
126       previous = sliceLocations[sliceIndex[i]];
127     }
128
129     //===========================================
130     // Create image  
131     gdcm::File * header = clitk::readDicomHeader(input_files[sliceIndex.front()]);
132     typedef itk::Image<float, 3> ImageType;
133     typedef itk::Image<signed short, 3> ImageTypeOutput;
134     ImageType::SpacingType spacing;
135     ImageType::SizeType    size;
136     ImageType::PointType    origin;
137     spacing[0] = header->GetXSpacing();
138     spacing[1] = header->GetYSpacing();
139     spacing[2] = dist;
140     size[0]    = header->GetXSize();
141     size[1]    = header->GetYSize();
142     size[2]    = sliceIndex.size();
143     ///Optional use of special origin scheme used at CLB
144     if (args_info.focal_origin_flag)
145     {
146         origin[0]  = -spacing[0]*size[0]/2;
147         origin[1]  = -spacing[1]*size[1]/2;
148     }
149     else
150     {
151         origin[0]  = header->GetXOrigin();
152         origin[1]  = header->GetYOrigin();
153     }
154     origin[2]  = header->GetZOrigin();
155     ImageTypeOutput::Pointer output = ImageTypeOutput::New();
156     itk::ImageRegion<3> region;
157     region.SetSize(size);
158     output->SetRegions(region);
159     output->SetOrigin(origin);
160     output->SetSpacing(spacing);
161     output->Allocate();
162   
163     //===========================================
164     // Fill image
165     ImageType::Pointer slice;
166     typedef itk::ImageRegionIterator<ImageTypeOutput> IteratorType;
167     typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
168     IteratorType po(output, output->GetLargestPossibleRegion());
169     po.Begin();
170     for(unsigned int i=0; i<sliceIndex.size(); i++) {
171       slice = clitk::readImage<ImageType>(input_files[sliceIndex[i]]);
172       ConstIteratorType pi(slice, slice->GetLargestPossibleRegion());
173       while ( !pi.IsAtEnd() ) {
174         po.Set((signed short)pi.Get());
175         ++po;
176         ++pi;
177       }
178     }
179
180     //===========================================
181     // Write image
182     clitk::writeImage<ImageTypeOutput>(output, 
183                                        args_info.output_arg, 
184                                        args_info.verbose_flag);
185   
186     // this is the end my friend 
187     return 0;
188 }