]> Creatis software - clitk.git/blob - clitkImageExtractLine.cxx
b46a5d06a97724d660f82be9d3ca387f1c50b5ef
[clitk.git] / clitkImageExtractLine.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 CLITKIMAGEEXTRACTLINE_CXX
19 #define CLITKIMAGEEXTRACTLINE_CXX
20 /**
21    -------------------------------------------------
22    * @file   clitkImageExtractLine.cxx
23    * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
24    * @date   23 Feb 2008 08:37:53
25    * @modified by Loïc Grevillot <Loic.Grevillot@creatis.insa-lyon.fr>
26    * @date   10 March 2011
27       * Option -I added, in order to integrate plans perpendicular to a line
28    
29    -------------------------------------------------*/
30
31 // clitk include
32 #include "clitkImageExtractLine_ggo.h"
33 #include "clitkIO.h"
34 #include "clitkImageCommon.h"
35 #include "clitkCommon.h"
36 #include <itkLineConstIterator.h>
37
38 //--------------------------------------------------------------------
39 int main(int argc, char * argv[])
40 {
41
42   // Init command line
43   GGO(clitkImageExtractLine, args_info);
44   CLITK_INIT;
45
46   // Declare main types
47   typedef float PixelType;
48   const unsigned int Dimension=3;
49   typedef itk::Image<PixelType, Dimension> ImageType;
50   typedef itk::Size<Dimension> SizeType;
51
52   // Check options
53   if (args_info.firstIndex_given != Dimension) {
54     std::cerr << "Please give " << Dimension << "values to --firstIndex option" << std::endl;
55     exit(0);
56   }
57   if (args_info.lastIndex_given != Dimension) {
58     std::cerr << "Please give " << Dimension << "values to --lastIndex option" << std::endl;
59     exit(0);
60   }
61
62   // Read image
63   ImageType::Pointer input = clitk::readImage<ImageType>(args_info.input_arg, args_info.verbose_flag);
64
65   // Get first and last index
66   typedef ImageType::IndexType IndexType;
67   IndexType firstIndex;
68   IndexType lastIndex;
69   ImageType::SpacingType spacing = input->GetSpacing();
70   double length = 0.0;
71   for(unsigned int i=0; i<Dimension; i++) {
72     firstIndex[i] = args_info.firstIndex_arg[i];
73     lastIndex[i] = args_info.lastIndex_arg[i];
74     if (args_info.mm_flag) {
75       firstIndex[i] /= spacing[i];
76       lastIndex[i] /= spacing[i];
77     }
78     length += pow(lastIndex[i]*spacing[i]-firstIndex[i]*spacing[i],2);
79   }
80   length = sqrt(length);
81
82   // Loop
83 //   std::vector<double> depth;
84 //   std::vector<double> values;
85 //   itk::LineConstIterator<ImageType> iter(input, firstIndex, lastIndex);
86 //   iter.GoToBegin();
87 //   while (!iter.IsAtEnd()) {
88 //     values.push_back(iter.Get());
89 //     ++iter;
90 //   }
91 //   double step = length/values.size();
92
93   std::vector<double> depth;
94   std::vector<double> values;
95   itk::LineConstIterator<ImageType> iter(input, firstIndex, lastIndex);
96   int direction=0;
97   
98   // args_info.integral_arg=0, so, it does not compute the integral
99   if (args_info.integral_arg==0){
100     iter.GoToBegin();
101     while (!iter.IsAtEnd()) {
102       values.push_back(iter.Get());
103       ++iter;
104     }
105   }
106   // args_info.integral_arg=1, so, it computes the integral
107   if (args_info.integral_arg!=0){
108     int a=0, b=0, c=0;
109     if (args_info.integralAxis_arg==0){
110         a=1;
111         b=0;
112         c=2;
113     }
114     else if (args_info.integralAxis_arg==1){
115         a=0;
116         b=1;
117         c=2;
118     }
119     else if (args_info.integralAxis_arg==2){
120         a=2;
121         b=0;
122         c=1;
123     }
124     else {std::cout<<"Wrong axis"<<std::endl;
125       exit(0);
126     }
127   
128   length=(lastIndex[b]-firstIndex[b])*spacing[b];
129   
130   std::cout<<"The line is extracted along axis "<<args_info.integralAxis_arg<<std::endl;
131   std::cout<<"The line is integrated between "<<args_info.firstIndex_arg[a]<<" and "<<args_info.lastIndex_arg[a]<<" along axis "<<a<<std::endl;
132   std::cout<<"The line is integrated between "<<args_info.firstIndex_arg[c]<<" and "<<args_info.lastIndex_arg[c]<<" along axis "<<c<<std::endl;
133   
134       SizeType dim;
135       dim=input->GetLargestPossibleRegion().GetSize();
136       DD(dim);
137       DD(direction);
138       
139       std::vector<double> val(dim[b]);
140         for (uint i=0; i<dim[b]; i++)
141           val[i]=0;
142         
143       int k;
144         
145       for (int i=args_info.firstIndex_arg[a]; i<args_info.lastIndex_arg[a]; i++){
146         for (int j=args_info.firstIndex_arg[c]; j<args_info.lastIndex_arg[c]; j++){
147     //      std::cout<<"i "<<i<<"  j "<<j<<std::endl;
148           k=0;
149           firstIndex[a]=i;
150           firstIndex[c]=j;
151           lastIndex[a]=i;
152           lastIndex[c]=j;
153     //      std::cout<<"A"<<std::endl;
154           itk::LineConstIterator<ImageType> iter(input, firstIndex, lastIndex);
155           iter.GoToBegin();
156     //      std::cout<<"B"<<std::endl;
157           val[k]+=iter.Get();
158           k++;
159     //      std::cout<<"C"<<std::endl;
160           while (!iter.IsAtEnd()) {
161     //  std::cout<<"D "<<k<<std::endl;
162             val[k]+=iter.Get();
163             ++iter;
164             k++;
165           }
166         }
167       }
168   
169       for (unsigned int i=0; i<dim[b]; i++){
170         values.push_back(val[i]);
171       }
172   }
173   
174   double step = length/values.size();
175   DD(values.size());
176   
177   // If isocenter is used
178   double isoDistance = 0.0;
179   if (args_info.isocenter_given) { // isoCenter is in mm
180     IndexType isoCenter;
181     for(unsigned int i=0; i<Dimension; i++) {
182       isoCenter[i] = args_info.isocenter_arg[i];
183       isoDistance += pow(isoCenter[i] - firstIndex[i]*spacing[i],2);
184     }
185     DD(isoCenter);
186     isoDistance = sqrt(isoDistance);
187     DD(isoDistance);
188   }
189
190   // Write result
191   std::ofstream os(args_info.output_arg);
192   os << "# clitkImageExtractLine from " << args_info.input_arg << std::endl
193      << "# \t firstIndex = " << firstIndex << std::endl
194      << "# \t lastIndex  = " << lastIndex << std::endl
195      << "# \t nb values  = " << values.size() << std::endl
196      << "# \t length   = " << length << std::endl
197      << "# \t step       = " << step << std::endl;
198   if (args_info.depth_flag) {
199     double lg = -isoDistance;
200     for(unsigned int i=0; i<values.size(); i++) {
201       os << lg << " " << values[i] << std::endl;
202       lg += step;
203     }
204     os.close();
205   } else {
206     for(unsigned int i=0; i<values.size(); i++) {
207       os << values[i] << std::endl;
208     }
209     os.close();
210   }
211
212   // this is the end my friend
213   return 0;
214 } // end main
215
216 #endif //define CLITKIMAGEEXTRACTLINE_CXX