tools/clitkLineProfileGenericFilter.txx

   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 clitkLineProfileGenericFilter_cxx
  19 #define clitkLineProfileGenericFilter_cxx
  20
  21 /* =================================================
  22  * @file   clitkLineProfileGenericFilter.cxx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "itkBresenhamLine.h"
  31
  32 namespace clitk
  33 {
  34
  35 //-----------------------------------------------------------
  36   // Update
  37   //-----------------------------------------------------------
  38   template<class InputImageType>
  39   void LineProfileGenericFilter::UpdateWithInputImageType()
  40   {
  41     typedef InputImageType ImageType;
  42
  43     if (m_Verbose)
  44       std::cout << "LineProfileGenericFilter::UpdateWithInputImageType" << std::endl;
  45
  46     // type checks
  47     if (m_ArgsInfo.p0_given != ImageType::ImageDimension ||
  48         m_ArgsInfo.p1_given != ImageType::ImageDimension)
  49         throw std::logic_error("Dimension of input points and image do not match.");
  50
  51     typename ImageType::Pointer image = this->template GetInput<InputImageType>(0);
  52     typename ImageType::RegionType region = image->GetLargestPossibleRegion();
  53
  54     typedef typename ImageType::PointType PointType;
  55     PointType p0, p1;
  56     for (unsigned int i = 0; i < ImageType::ImageDimension; i++) {
  57       p0[i] = m_ArgsInfo.p0_arg[i];
  58       p1[i] = m_ArgsInfo.p1_arg[i];
  59     }
  60
  61     // compute params of line between p0 and p1
  62     // length (magnitude) must be an integer value, so it's
  63     // the max distance along one the axes plus one to account
  64     // for the last point.
  65     typedef itk::BresenhamLine<ImageType::ImageDimension> LineType;
  66     typename LineType::LType direction = p1 - p0;
  67     typename LineType::LType::RealValueType mag = 0;
  68     for (unsigned int i = 0; i < ImageType::ImageDimension; i++) {
  69       if (direction[i] > mag)
  70         mag = direction[i];
  71     }
  72     mag++;
  73
  74     if (m_Verbose)
  75       std::cout << "Building line with direction = " << direction << " and length = " << mag << "..." << std::endl;
  76
  77     // build the line itself
  78     LineType line;
  79     typename LineType::OffsetArray offsets;
  80     offsets = line.BuildLine(direction, mag);
  81
  82     if (m_Verbose)
  83       std::cout << "Line has " << offsets.size() << " points" << std::endl;
  84
  85     // fill the output vectors
  86     typedef typename ImageType::PixelType PixelType;
  87     typedef typename ImageType::OffsetType OffsetType;
  88     typedef typename ImageType::IndexType IndexType;
  89     typedef std::vector<IndexType> IndexArrayType;
  90     typedef std::vector<PixelType> ValueArrayType;
  91
  92     IndexType index0, index1;
  93     for (unsigned int i = 0; i < ImageType::ImageDimension; i++) {
  94       index0[i] = m_ArgsInfo.p0_arg[i];
  95       index1[i] = m_ArgsInfo.p1_arg[i];
  96     }
  97
  98     if (m_Verbose)
  99       std::cout << "Getting profile along line..." << std::endl;
 100
 101     IndexType index;
 102     IndexArrayType line_indices;
 103     ValueArrayType line_values;
 104     for (size_t i = 0; i < offsets.size(); i++)
 105     {
 106       index = index0 + offsets[i];
 107       if (m_Verbose) {
 108         std::cout << "index " << i << " = " << index << std::endl;
 109       }
 110
 111       if (region.IsInside(index)) {
 112         if (m_Verbose)
 113           std::cout << "value " << i << " = " << image->GetPixel(index) << std::endl;
 114
 115         line_indices.push_back(index);
 116         line_values.push_back(image->GetPixel(index));
 117       }
 118       else if (m_Verbose)
 119         std::cout << "index outside image" << std::endl;
 120     }
 121
 122     if (m_Verbose) {
 123       std::cout << "I bring to you The Computed Points!" << std::endl;
 124     }
 125
 126     for (size_t i = 0; i < line_indices.size(); i++) {
 127       std::cout << line_indices[i] << " " << line_values[i] << std::endl;
 128     }
 129   }
 130
 131
 132 } //end clitk
 133
 134 #endif  //#define clitkLineProfileGenericFilter_cxx