]> Creatis software - clitk.git/blob - tools/clitkImageExtractLine.cxx
Romulo:
[clitk.git] / tools / 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://oncora1.lyon.fnclcc.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   else if (lastIndex[0]-firstIndex[0]==0 && lastIndex[1]-firstIndex[1]==0 && lastIndex[2]-firstIndex[2]>0)
108     direction=1;
109   else if (lastIndex[0]-firstIndex[0]==0 && lastIndex[1]-firstIndex[1]>0 && lastIndex[2]-firstIndex[2]==0)
110     direction=2;
111   else if (lastIndex[0]-firstIndex[0]>0 && lastIndex[1]-firstIndex[1]==0 && lastIndex[2]-firstIndex[2]==0)
112     direction=3;
113   else{
114     //std::cout<<lastIndex[0]-firstIndex[0]<<"  "<<lastIndex[1]-firstIndex[1]<<"  "<<lastIndex[3]-firstIndex[3]<<std::endl;
115     std::cout<<"Index are not defined along a straight along x or y or z axis."<<std::endl;
116     std::cout<<"The line cannot be extracted."<<std::endl;
117     exit(0);
118   }
119   
120   if (args_info.integral_arg!=0){
121     SizeType dim;
122     dim=input->GetLargestPossibleRegion().GetSize();
123     DD(dim);
124     DD(direction);
125     
126     int a=0, b=0, c=0;
127     
128     if (direction==2){
129       a=0;
130       b=1;
131       c=2;
132     }
133     if (direction==1){
134       a=0;
135       b=2;
136       c=1;
137     }    
138     if (direction==3){
139       a=2;
140       b=0;
141       c=1;
142     } 
143     
144     double val[dim[b]];
145       for (int i=0; i<dim[b]; i++)
146         val[i]=0;
147       
148     int k;
149     for (int i=0; i<dim[a]; i++){
150       for (int j=0; j<dim[c]; j++){
151   //      std::cout<<"i "<<i<<"  j "<<j<<std::endl;
152         k=0;
153         firstIndex[a]=i;
154         firstIndex[c]=j;
155         lastIndex[a]=i;
156         lastIndex[c]=j;
157   //      std::cout<<"A"<<std::endl;
158         itk::LineConstIterator<ImageType> iter(input, firstIndex, lastIndex);
159         iter.GoToBegin();
160   //      std::cout<<"B"<<std::endl;
161         val[k]+=iter.Get();
162         k++;
163   //      std::cout<<"C"<<std::endl;
164         while (!iter.IsAtEnd()) {
165   //    std::cout<<"D "<<k<<std::endl;
166           val[k]+=iter.Get();
167           ++iter;
168           k++;
169         }
170       }
171     }
172   
173     for (unsigned int i=0; i<dim[b]; i++){
174       values.push_back(val[i]);
175     }
176   }
177   
178   double step = length/values.size();
179   DD(values.size());
180   
181   // If isocenter is used
182   double isoDistance = 0.0;
183   if (args_info.isocenter_given) { // isoCenter is in mm
184     IndexType isoCenter;
185     for(unsigned int i=0; i<Dimension; i++) {
186       isoCenter[i] = args_info.isocenter_arg[i];
187       isoDistance += pow(isoCenter[i] - firstIndex[i]*spacing[i],2);
188     }
189     DD(isoCenter);
190     isoDistance = sqrt(isoDistance);
191     DD(isoDistance);
192   }
193
194   // Write result
195   std::ofstream os(args_info.output_arg);
196   os << "# clitkImageExtractLine from " << args_info.input_arg << std::endl
197      << "# \t firstIndex = " << firstIndex << std::endl
198      << "# \t lastIndex  = " << lastIndex << std::endl
199      << "# \t nb values  = " << values.size() << std::endl
200      << "# \t length   = " << length << std::endl
201      << "# \t step       = " << step << std::endl;
202   if (args_info.depth_flag) {
203     double lg = -isoDistance;
204     for(unsigned int i=0; i<values.size(); i++) {
205       os << lg << " " << values[i] << std::endl;
206       lg += step;
207     }
208     os.close();
209   } else {
210     for(unsigned int i=0; i<values.size(); i++) {
211       os << values[i] << std::endl;
212     }
213     os.close();
214   }
215
216   // this is the end my friend
217   return 0;
218 } // end main
219
220 #endif //define CLITKIMAGEEXTRACTLINE_CXX