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