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