tools/clitkProfileImageGenericFilter.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 clitkProfileImageGenericFilter_cxx
  19 #define clitkProfileImageGenericFilter_cxx
  20
  21 /* =================================================
  22  * @file   clitkProfileImageGenericFilter.cxx
  23  * @author Thomas Baudier <thomas.baudier@creatis.insa-lyon.fr>
  24  * @date   22 dec 2015
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "clitkProfileImageGenericFilter.h"
  31
  32 // itk include
  33 #include <itkLineIterator.h>
  34 #include <itkPoint.h>
  35
  36 #include <clitkCommon.h>
  37
  38
  39
  40 namespace clitk
  41 {
  42
  43 //--------------------------------------------------------------------
  44 ProfileImageGenericFilter::ProfileImageGenericFilter():
  45   ImageToImageGenericFilter<Self>("ProfileImage")
  46 {
  47   InitializeImageType<2>();
  48   InitializeImageType<3>();
  49   InitializeImageType<4>();
  50 }
  51 //--------------------------------------------------------------------
  52
  53
  54 //--------------------------------------------------------------------
  55 template<unsigned int Dim>
  56 void ProfileImageGenericFilter::InitializeImageType()
  57 {
  58   ADD_DEFAULT_IMAGE_TYPES(Dim);
  59 }
  60 //--------------------------------------------------------------------
  61
  62
  63 //--------------------------------------------------------------------
  64 vtkFloatArray* ProfileImageGenericFilter::GetArrayX()
  65 {
  66   return(mArrayX);
  67 }
  68 //--------------------------------------------------------------------
  69
  70
  71 //--------------------------------------------------------------------
  72 vtkFloatArray* ProfileImageGenericFilter::GetArrayY()
  73 {
  74   return(mArrayY);
  75 }
  76 //--------------------------------------------------------------------
  77
  78
  79 //--------------------------------------------------------------------
  80 vtkFloatArray* ProfileImageGenericFilter::GetCoord()
  81 {
  82   return(mCoord);
  83 }
  84 //--------------------------------------------------------------------
  85
  86
  87 //--------------------------------------------------------------------
  88 void ProfileImageGenericFilter::SetArgsInfo(const args_info_type & a)
  89 {
  90   mArgsInfo=a;
  91   if (mArgsInfo.verbose_given)
  92     SetIOVerbose(mArgsInfo.verbose_flag);
  93   if (mArgsInfo.imagetypes_given && mArgsInfo.imagetypes_flag)
  94     this->PrintAvailableImageTypes();
  95
  96   if (mArgsInfo.input_given) {
  97     SetInputFilename(mArgsInfo.input_arg);
  98   }
  99   if (mArgsInfo.output_given) {
 100     SetOutputFilename(mArgsInfo.output_arg);
 101   }
 102 }
 103 //--------------------------------------------------------------------
 104
 105
 106 //--------------------------------------------------------------------
 107 // Update with the number of dimensions and the pixeltype
 108 //--------------------------------------------------------------------
 109 template<class InputImageType>
 110 void
 111 ProfileImageGenericFilter::UpdateWithInputImageType()
 112 {
 113
 114   // Reading input
 115   typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
 116   typedef typename InputImageType::PixelType PixelType;
 117   typedef typename InputImageType::IndexType IndexType;
 118
 119   mArrayX = vtkSmartPointer<vtkFloatArray>::New();
 120   mArrayY = vtkSmartPointer<vtkFloatArray>::New();
 121   mCoord = vtkSmartPointer<vtkFloatArray>::New();
 122   mCoord->SetNumberOfComponents(InputImageType::ImageDimension);
 123   mCoordmm = vtkSmartPointer<vtkFloatArray>::New();
 124   mCoordmm->SetNumberOfComponents(InputImageType::ImageDimension);
 125   mDimension = InputImageType::ImageDimension;
 126
 127   /*typename InputImageType::Pointer outputImage;
 128   outputImage = InputImageType::New();
 129
 130   outputImage->SetRegions(input->GetLargestPossibleRegion());
 131   outputImage->Allocate();
 132   outputImage->FillBuffer(0); */
 133
 134   //Iterator
 135   IndexType pointBegin, pointEnd;
 136
 137   for (int i = 0; i < mArgsInfo.point1_given; ++i) {
 138     pointBegin[i] = mArgsInfo.point1_arg[i];
 139     pointEnd[i] = mArgsInfo.point2_arg[i];
 140   }
 141
 142   itk::LineConstIterator<InputImageType> itProfile(input, pointBegin, pointEnd);
 143   itProfile.GoToBegin();
 144   int lineNumber(1);
 145   double *tuple;
 146   double distance;
 147   tuple = new double[InputImageType::ImageDimension];
 148   itk::Point<double, InputImageType::ImageDimension> transformedFirstPoint;
 149   itk::Point<double, InputImageType::ImageDimension> transformedCurrentPoint;
 150
 151   input->TransformIndexToPhysicalPoint(itProfile.GetIndex(), transformedFirstPoint);
 152
 153   while (!itProfile.IsAtEnd())
 154   {
 155     // Fill in the table the intensity value
 156     mArrayY->InsertNextTuple1(itProfile.Get());
 157
 158     for (int i=0; i<InputImageType::ImageDimension; ++i) {
 159         tuple[i] = itProfile.GetIndex()[i];
 160     }
 161
 162     input->TransformIndexToPhysicalPoint(itProfile.GetIndex(), transformedCurrentPoint);
 163     distance = transformedFirstPoint.EuclideanDistanceTo(transformedCurrentPoint);
 164
 165     // Fill in the table the distance value
 166     mArrayX->InsertNextTuple1(distance);
 167
 168     // Fill in the table the voxel coordinate value
 169     mCoord->InsertNextTuple(tuple); //index
 170     for (int i=0; i<InputImageType::ImageDimension; ++i) {
 171         tuple[i] = transformedCurrentPoint[i];
 172     }
 173     mCoordmm->InsertNextTuple(tuple); //mm
 174     ++lineNumber;
 175     ++itProfile;
 176   }
 177
 178   if (mArgsInfo.output_given) {
 179     std::string str(mArgsInfo.output_arg);
 180     this->WriteOutput(str);
 181   }
 182
 183   /*
 184   itk::LineIterator<InputImageType> otProfile(outputImage, pointBegin, pointEnd);
 185   otProfile.GoToBegin();
 186   while (!otProfile.IsAtEnd())
 187   {
 188     otProfile.Set(1.0);
 189     ++otProfile;
 190   }
 191
 192   this->template SetNextOutput<InputImageType>(outputImage): */
 193
 194   delete [] tuple;
 195 }
 196 //--------------------------------------------------------------------
 197
 198
 199 //--------------------------------------------------------------------
 200 void ProfileImageGenericFilter::WriteOutput(std::string outputFilename)
 201 {
 202   ofstream fileOpen(outputFilename.c_str(), std::ofstream::trunc);
 203
 204   if(!fileOpen) {
 205       cerr << "Error during saving" << endl;
 206       return;
 207   }
 208
 209   double *tuple;
 210   tuple = new double[mDimension];
 211   int i(0);
 212   fileOpen << "The Bresenham algorithm is used to travel along the line. Values represent the center of each crossed voxel (in voxel and mm)" << endl;
 213   fileOpen << "Id" << "\t" << "Value" << "\t" ;
 214   fileOpen << "x(vox)" << "\t" << "y(vox)" << "\t";
 215   if (mDimension >=3)
 216       fileOpen << "z(vox)" << "\t";
 217   if (mDimension >=4)
 218       fileOpen << "t" << "\t";
 219   fileOpen << "x(mm)" << "\t" << "y(mm)" << "\t";
 220   if (mDimension >=3)
 221       fileOpen << "z(mm)" << "\t";
 222   if (mDimension >=4)
 223       fileOpen << "t" << "\t";
 224   fileOpen << endl;
 225
 226   while (i<mArrayX->GetNumberOfTuples()) {
 227       fileOpen << i << "\t" << mArrayY->GetTuple(i)[0] << "\t" ;
 228
 229       mCoord->GetTuple(i, tuple);
 230       for (int j=0; j<mDimension ; ++j) {
 231           fileOpen << tuple[j] << "\t" ;
 232       }
 233       mCoordmm->GetTuple(i, tuple);
 234       for (int j=0; j<mDimension ; ++j) {
 235           fileOpen << tuple[j] << "\t" ;
 236       }
 237       if (mDimension == 4) {
 238           fileOpen << tuple[3] << "\t" ;
 239       }
 240       fileOpen << endl;
 241       ++i;
 242   }
 243
 244   delete [] tuple;
 245
 246   fileOpen.close();
 247 }
 248 //--------------------------------------------------------------------
 249
 250
 251 }//end clitk
 252
 253 #endif  //#define clitkProfileImageGenericFilter_cxx